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V.  Introduction 


An  investigation  of  three-dimensional  ultrasonic  mammography  is  underway.  The  goal 
of  the  research  is  improved  diagnosis  of  breast  cancer  by  quantitative,  high-resolution  three- 
dimensional  ultrasonic  imaging.  This  goal  is  being  reached  by  a  thorough  program  that 
synthesizes  recent  advances  in  tissue  modeling,  adaptive  imaging,  instrumentation,  and  sig¬ 
nal  processing.  The  final  result  of  the  research  will  be  a  major  advance  in  quantitative 
three-dimensional  ultrasonic  mammography.  Improved  resolution,  accurate  quantitative  in¬ 
formation  on  tissue  properties,  and  precise  determination  of  three-dimensional  breast  struc¬ 
ture  will  provide  crucial  new  information  for  detection,  diagnosis,  and  monitoring  of  breast 
cancer.  The  goal  of  three-dimensional  quantitative  imaging  is  currently  being  achieved  using 
novel  inverse  scattering  methods  invented  by  the  Principal  Investigator  and  coworkers.  Use  of 
full  time-domain  scattering  information  provides  images  with  high  point  resolution,  contrast 
resolution,  and  quantitative  accuracy  without  significant  artifacts.  Nonlinear  forms  of  these 
methods  provide  a  robust  approach  to  adaptive  imaging  that  is  based  on  compensation  for 
three-dimensional  scattering  from  actual  tissue  structure.  Unlike  previous  adaptive  imaging 
methods  based  on  assumptions  of  phase-screen  aberrators  and  point  scatterers,  these  meth¬ 
ods  provide  aberration  correction  ideally  suited  to  distributed  inhomogeneous  tissue  like  the 
breast.  A  unique  and  innovative  aspect  of  the  research  is  the  use  of  realistic  tissue  models 
for  ultrasonic  propagation  through  breast  tissue.  Such  models  have  been  shown  to  realisti¬ 
cally  model  wavefront  distortion  in  the  human  abdominal  wall,  but  have  not  to  date  been 
applied  to  the  human  breast.  Tissue  modeling  techniques  employ  tissue  maps  obtained  from 
specimen  cross  sections  as  well  as  from  newly  available  high-resolution  volume  photographic 
data.  Calculated  scattering  from  these  tissue  models  will  provide  accurate  characterization 
of  ultrasonic  propagation  within  breast  tissue  and  will  also  provide  realistic  data  for  quanti¬ 
tative  imaging  algorithms.  The  above  studies  will  facilitate  the  application  of  breakthroughs 
from  tissue  modeling,  inverse  scattering,  and  signal  processing  to  the  critical  application 
of  ultrasonic  mammography.  Completion  of  the  proposed  research  will  make  possible  new 
mammographic  applications  of  ultrasound  that  will  provide  clinicians  with  previously  un¬ 
available  information  and  detail.  The  end  result  will  be  a  lower-cost,  more  effective,  and 
safer  modality  for  diagnosis,  detection,  and  monitoring  of  breast  cancer. 
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VI.  Body  of  Report 

Below,  the  accomplishments  of  the  first  year  for  this  project  are  summarized  under  the 
four  categories  of  the  approved  Statement  of  Work.  Details  are  given  for  how  accomplish¬ 
ments  to  date  fit  into  the  overall  research  plan.  Where  applicable,  brief  descriptions  of 
planned  research  indicate  how  the  remainder  of  the  Statement  of  Work  will  be  fulfilled. 

A.  Quantitative  Imaging  Algorithm  Development 

In  the  first  year  of  this  USAMRMC-funded  Breast  Cancer  Research  Program  project, 
great  progress  has  been  made  on  development  of  new  quantitative  imaging  algorithms  for 
three-dimensional  ultrasonic  mammography.  Research  on  algorithm  development  has  oc¬ 
curred  in  two  areas:  (1)  a  new  time-domain  diffraction  tomography  method  for  wideband 
quantitative  imaging  and  (2)  nonlinear  extensions  to  a  novel  eigenfunction-based  inverse 
scattering  method. 

The  new  time-domain  diffraction  tomography  method  [1,2],  which  is  specifically  designed 
for  ultrasonic  mammography,  allows  quantitative,  high- resolution  images  to  be  obtained  us¬ 
ing  direct  synthetic-aperture  processing  of  time-domain  scattered  fields.  This  method  pro¬ 
vides  tomographic  images  of  inhomogeneous  media  using  scattering  measurements  made  on 
a  surface  surrounding  the  medium  of  interest,  e.g.,  on  a  circle  for  two-dimensional  problems 
or  on  a  sphere  for  three-dimensional  problems.  Images  of  compressibility  variations  are  then 
reconstructed  using  coherent  combination  of  the  far-field  scattered  waveforms,  delayed  and 
summed  in  a  manner  that  numerically  focuses  on  the  unknown  medium.  This  approach 
is  closely  related  to  synthetic  aperture  imaging;  however,  unlike  conventional  synthetic- 
aperture  methods,  the  present  method  provides  quantitative  reconstructions  of  compressibil¬ 
ity  variations,  analogous  to  frequency-compounded  filtered  backpropagation  images  weighted 
by  the  spectrum  of  the  incident  wave.  A  complete  description  of  the  method,  with  example 
two-dimensional  and  three-dimensional  reconstruction  results,  is  given  in  Ref.  [2],  included 
as  Appendix  A  of  this  report.  A  brief  abstract  describing  the  work,  presented  at  Forum 
Acusticum  99,  is  included  in  Appendix  D. 

The  results  of  the  new  time-domain  imaging  method  are  exciting  for  several  reasons. 
First,  the  images  are  both  higher  in  quality  and  more  efficiently  computed  than  conventional 
single-frequency  quantitative  images.  The  high  point  and  contrast  resolution,  as  well  as 
absence  of  artifacts  usually  associated  with  diffraction  tomography,  suggests  that  this  method 
will  be  very  useful  for  detection  and  characterization  of  breast  lesions.  Second,  because  of 
the  close  analogy  between  the  new  method  and  delay-and-sum  imaging,  the  new  method 
could  be  implemented  in  hardware  using  beamforming  technology  already  present  on  digital 
ultrasound  scanners.  Furthermore,  the  new  time-domain  method  can  easily  incorporate 
other  imaging  techniques  (e.g.,  time-gain  compensation,  harmonic  imaging,  and  aberration 
correction)  currently  used  in  clinical  and  experimental  B-scan  systems. 

Work  was  also  performed  on  nonlinear  (aberration-corrected)  implementation  of  a 
method  for  quantitative  imaging  using  eigenfunctions  of  the  scattering  operator  [3]-[4].  A 
1998  presentation  of  this  method  at  a  meeting  of  the  American  Institute  of  Ultrasound  in 
Medicine  [5]  was  received  with  interest.  Preliminary  work  on  iterative  nonlinear  inverse 
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scattering,  with  application  to  quantitative  imaging  of  strongly  scattering  media  such  as  the 
human  breast,  shows  promise  for  the  potential  of  aberration  correction  to  improve  quantita¬ 
tive  ultrasonic  mammography. 

Work  on  algorithm  development  will  continue  with  synthesis  of  wideband  time-domain 
inverse  scattering  together  with  aberration  correction  based  on  quantitative  estimates  of 
the  unknown  medium.  Since  these  estimates  are  obtained  directly  from  inverse  scattering 
reconstructions,  iterative  reconstruction  procedures  will  provide  the  highest  quantitative 
accuracy  available  for  the  data.  Since  the  nonlinear  eigenfunction  method  provides  a  simple 
approach  to  aberration  correction  (improved  reconstructions  are  obtained  using  numerical 
retransmission  of  eigenfunctions  into  the  estimated  medium),  this  procedure  is  much  more 
efficient  than  other  available  nonlinear  inverse  scattering  methods.  This  form  of  nonlinear 
aberration  correction,  in  conjunction  with  the  new  time-domain  imaging  method  described 
above,  provide  all  the  capabilities  required  for  the  proposed  method  of  three-dimensional 
quantitative  ultrasonic  mammography. 

B.  Tissue  Modeling 

Progress  toward  improved  scattering  models  for  ultrasound-breast  tissue  interaction  has 
been  made  in  several  related  studies. 

A  finite-difference  time-domain  simulation  of  ultrasonic  propagation  through  cross- 
sectional  models  of  chest  wall  tissue  [6,  7]  has  generated  several  important  results  that  are 
directly  relevant  to  the  planned  simulation  of  propagation  through  breast  tissue.  First,  this 
work  expands  upon  previous  models  [8,  9]  by  including  tissue-dependent  absorption  effects. 
The  absorption  model  developed  will  be  an  important  part  of  the  planned  two-dimensional 
and  three-dimensional  simulations  of  ultrasonic  scattering  from  breast  tissue.  Second,  the 
chest  study  contained  an  analysis  of  the  frequency  dependence  of  ultrasonic  wavefront  dis¬ 
tortion.  A  comparison  of  distortion  effects  for  varying  pulse  center  frequencies  showed  that, 
for  soft  tissue  paths  through  the  chest  wall,  energy  level  and  waveform  distortion  increase 
markedly  with  rising  ultrasonic  frequency  and  that  arrival-time  fluctuations  increase  to  a 
lesser  degree.  Since  breast  tissue  has  been  observed  to  be  more  strongly  scattering  than 
other  soft  tissues  [10],  these  results  are  of  particular  importance  for  ultrasonic  mammogra¬ 
phy.  A  full  description  of  this  study,  with  numerical  results  including  statistical  summaries, 
is  given  in  Ref.  [7],  included  here  as  Appendix  B.  An  abstract  describing  the  work,  presented 
at  the  136th  Meeting  of  the  Acoustical  Society  of  America,  is  included  in  Appendix  D. 

In  accordance  with  the  Statement  of  Work,  scattering  has  been  computed  for  canonically- 
shaped  inhomogeneities  including  spheres,  cylinders,  and  rectangular  slabs.  Time-domain 
scattered  fields  were  computed  using  exact  solutions  [12]  for  spheres  and  cylinders  with  sound 
speed  and  density  different  from  the  surrounding  medium.  Time-domain  scattered  fields  were 
also  computed,  using  a  weak  scattering  approximation,  for  rectangular  slabs  with  a  sound 
speed  different  from  the  surrounding  medium.  The  purpose  of  these  computations  was  to 
provide  definitive  test  data  for  quantitative  imaging  algorithms.  Example  reconstructions 
based  on  these  computations  are  provided  in  Ref.  [2],  provided  as  Appendix  A. 

A  new  fc-space  method  was  implemented  for  simulating  wave  propagation  in  soft  tissue 
[11].  Numerical  results,  including  quantitative  comparisons  of  accuracy,  showed  that  the  k- 
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space  method  can  achieve  very  accurate  results  for  large-scale  ultrasonic  propagation  through 
soft  tissue.  The  computational  efficiency  obtained  with  the  new  k- space  method  is  much 
greater  than  that  of  existing  finite  difference  or  pseudospectral  methods.  The  method  is 
described  and  quantitatively  analyzed  in  Ref.  [11],  included  here  as  Appendix  C.  Because 
of  the  remarkable  accuracy  and  efficiency  of  the  new  /c-space  method,  a  three-dimensional 
version  of  this  method  will  be  employed  for  three-dimensional  scattering  computations  in  the 
remainder  of  this  project.  A  three-dimensional  version  of  the  /c-space  code,  incorporating 
the  absorption  model  of  Ref.  [7]  (Appendix  B)  as  well  as  absorbing  boundary  conditions, 
has  been  implemented  and  is  ready  for  use  with  breast  tissue  models.  The  existing  code  is 
also  fully  parallelizable,  since  most  of  the  computational  load  is  associated  with  execution 
of  multithreaded  fast  Fourier  transform  computations. 

Work  is  currently  ongoing  for  optimization  of  speed  and  accuracy  of  the  /c-space  method 
for  scattering  computations.  Preliminary  comparisons  between  the  k- space  code  and  a  com¬ 
mercial  pseudospectral  code  have  shown  that  the  /c-space  code  allows  accurate  results  to  be 
obtained  with  much  larger  time  steps  and  fewer  independent  variables,  so  that  the  k- space 
code  is  superior  both  in  computational  speed  and  in  storage  requirements.  Both  of  these 
considerations  are  critical  for  the  planned  three-dimensional  computations  of  scattering  from 
breast  tissue. 

Now  that  a  robust,  efficient  three-dimensional  method  for  computation  of  ultrasonic  prop¬ 
agation  is  available,  work  on  mapping  of  breast  tissue  has  begun  in  earnest.  In  collaboration 
with  colleagues  at  the  University  of  Rochester,  cross-sectional  breast  tissue  specimens  have 
been  sectioned  and  stained  for  high-resolution  segmentation  by  tissue  type.  A  license  for 
use  of  the  Visible  Woman  data  set  has  been  obtained;  the  second  year  of  the  project  will 
include  a  substantial  effort  to  construct  three-dimensional  tissue  maps  from  the  provided 
photographic,  x-ray,  and  MRI  data. 

C.  Quantitative  Imaging  Algorithm  Implementation 

Implementation  of  quantitative  imaging  algorithms,  including  the  time-domain  diffrac¬ 
tion  tomography  method  and  the  nonlinear  eigenfunction  method  described  above,  has  to 
date  been  primarily  performed  using  simulated  data.  Quantitative  reconstructions  performed 
using  the  time-domain  diffraction  tomography  method  are  shown  in  Ref.  [2]  (Appendix  A) 
for  simulated  two-dimensional  and  three-dimensional  scattering  data.  These  results  are 
promising  for  ultrasonic  mammography.  As  discussed  below  in  section  D  (Evaluation  and 
Comparison  of  Results),  time-domain  quantitative  images  show  parametric  accuracy,  high 
resolution,  and  few  artifacts.  Effects  of  limited  scattering  data,  shown  in  Fig.  3  of  Appendix 
A  for  the  2D  case,  indicate  that  accurate  images  can  be  obtained  without  the  necessity  of 
apertures  that  entirely  enclose  the  breast. 

Experimental  testing  of  the  new  time-domain  quantitative  imaging  method,  in  collabo¬ 
ration  with  colleagues  at  the  University  of  Rochester,  is  now  beginning.  Tissue-mimicking 
phantoms  composed  of  agar  gel  with  glass  spheres  have  been  constructed.  An  anthropomor¬ 
phic  breast-mimicking  phantom,  developed  in  collaboration  between  General  Electric,  Uni¬ 
versity  of  Wisconsin,  and  University  of  Rochester,  is  also  available  for  testing.  Time-domain 
scattering  data,  to  be  taken  with  the  2048-element  University  of  Rochester  ring  transducer 
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[13]  will  be  used  as  input  for  testing  of  the  new  time-domain  quantitative  imaging  method. 
Both  2D  and  3D  reconstructions  will  be  performed. 

As  proposed  in  the  original  Statement  of  Work,  a  nonlinear,  aberration-corrected  version 
of  the  new  time-domain  inverse  scattering  method  will  also  be  implemented  and  tested  with 
simulated  data  (from  the  breast  tissue  models  described  above  in  section  B,  Tissue  Modeling) 
as  well  as  experimental  data.  The  same  data  will  be  used  as  input  for  standard  synthetic- 
aperture  imaging,  so  that  image  quality  and  accuracy  can  be  directly  compared. 

D.  Evaluation  and  Comparison  of  Results 

Analysis  of  the  time-domain  imaging  results  for  simulated  data  is  given  in  Ref.  [2]  (Ap¬ 
pendix  A).  These  results  indicate  that  high  quantitative  accuracy  can  be  achieved.  Compu¬ 
tations  of  point-spread  functions  (Fig.  2  in  Appendix  A)  show  that  the  time-domain  method 
yields  higher  point  and  contrast  resolution  than  single-frequency  diffraction  tomography.  For 
the  3D  case,  the  level  of  the  first  sidelobe  is  reduced  by  13  dB,  while  the  second  sidelobe  is 
reduced  by  18  dB.  Because  of  the  broadband  scattering  information  employed,  the  width  of 
the  main  lobe  indicates  point  resolution  of  features  smaller  than  one-half  wavelength  at  the 
center  frequency.  Evaluation  of  quantitative  accuracy  (Fig.  4)  shows  that  the  time-domain 
diffraction  tomography  method  provides  parametric  accuracy  similar  to  established  single¬ 
frequency  methods  while  showing  much  more  immunity  from  artifacts.  Similar  quantitative 
analysis  will  be  performed  for  images  obtained  using  the  nonlinear  aberration-corrected  in¬ 
verse  scattering  method  mentioned  above. 

Two-dimensional  and  three-dimensional  reconstructions  performed  using  the  new  time- 
domain  quantitative  imaging  method  will  be  directly  compared  to  analogous  images  from 
other  modalities  such  as  x-ray  computed  tomography  and  conventional  B-scan  ultrasound. 
Consultants  Dr.  Jon  Meilstrup,  Dr.  Claudia  Kasales,  and  Dr.  Carroll  Osgood  will  help  with 
this  portion  of  evaluation  and  comparison  of  results.  Although  the  original  statement  of 
work  defers  this  image  evaluation  task  until  the  third  year  of  the  project,  this  work  may 
begin  during  the  second  year  if  no  unanticipated  experimental  difficulties  are  encountered. 

Regarding  the  task  of  publishing  results  in  archival  journals,  the  first  year  of  the  project 
has  already  resulted  in  three  submitted  manuscripts.  Further  publications  planned  to  be 
prepared  during  the  next  year  will  be  based  on  experimental  testing  of  the  new  time-domain 
quantitative  imaging  method,  optimization  and  benchmarking  of  the  k- space  algorithm,  and 
scattering  computations  employing  realistic  breast  tissue  models. 
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VII.  Key  Research  Accomplishments 

The  key  research  accomplishments  to  date  in  this  project  can  be  summarized  as  follows: 

•  Development  of  a  new  time-domain  quantitative  imaging  method  designed  specifically 
for  ultrasonic  mammography. 

•  Numerical  implementation  and  testing  of  the  new  time-domain  imaging  method,  show¬ 
ing  that  this  method  provides  high  accuracy  with  greater  efficiency  than  previous  in¬ 
verse  scattering  methods. 

•  Implementation  of  a  nonlinear  aberration-corrected  eigenfunction  inverse  scattering 
method. 

•  Implementation  of  a  tissue-dependent  absorption  model  for  simulations  of  scattering 
and  propagation. 

•  Characterization  of  the  frequency  dependence  of  ultrasonic  scattering  from  human  soft 
tissues. 

•  Exact  computation  of  time-domain  scattering  from  simple  objects  for  testing  of  quan¬ 
titative  imaging  methods. 

•  Implementation  and  testing  of  a  new  /c-space  method  for  computation  of  scattering, 
indicating  that  the  method  is  accurate  and  extremely  efficient,  and  therefore  ideal  for 
the  proposed  3D  computations  of  scattering  from  breast  tissue. 

•  Extension  of  the  new  fc-space  method  to  include  tissue-dependent  absorption,  absorbing 
boundary  layers,  and  three-dimensional  scattering. 

•  Use  of  the  /c-space  method  to  compute  time-domain  scattering  data  for  testing  of 
quantitative  imaging  algorithms. 

•  Quantitative  analysis  of  inverse  scattering  results,  indicating  that  time-domain  recon¬ 
structions  provide  much  higher  point  resolution,  contrast  resolution,  and  freedom  from 
artifacts  than  single-frequency  reconstructions. 


VIII.  Reportable  Outcomes 

Reportable  outcomes  for  the  research  have,  to  date,  included  three  papers  submitted 
to  archival  journals  (Refs.  [2],  [7],  and  [11]),  as  well  as  two  published  abstracts  presented 
at  scientific  meetings  (Refs.  [1]  and  [6]).  All  of  these  publications,  as  well  as  a  current 
Curriculum  Vitae  for  the  Principal  Investigator,  are  included  below  in  the  Appendices. 
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IX.  Conclusions 


The  first  year  of  this  USAMRMC-funded  project  has  yielded  considerable  progress  toward 
the  goal  of  improved  diagnosis  of  breast  cancer  by  three-dimensional  ultrasonic  imaging. 
Several  breakthroughs  have  been  made  which  will  provide  a  solid  foundation  for  the  continued 
tissue  modeling  and  quantitative  imaging  research  planned  for  the  last  two  years  of  the 
project. 

A  major  breakthrough  has  been  the  invention  and  implementation  of  a  new  time-domain 
diffraction  tomography  method  designed  for  ultrasonic  mammography.  This  method  is  po¬ 
tentially  of  very  great  importance  for  breast  cancer  diagnosis  for  several  reasons:  (1)  images 
have  higher  quality  than  that  achievable  by  conventional  inverse  scattering  methods  or  by 
current  ultrasound  scanners,  (2)  tissue  parameters  are  computed  and  quantitatively  imaged 
with  high  accuracy,  and  (3)  the  close  analogy  between  the  new  method  and  conventional 
synthetic-aperture  imaging  will  allow  rapid  implementation  of  the  new  method  on  hardware 
similar  to  currently  used  beamformers.  The  additional  improvement  of  aberration  correction 
should  increase  the  value  of  this  method  even  more,  because  the  strong  scattering  inherent 
to  breast  tissue  is  an  important  limiting  factor  to  existing  ultrasonic  imaging  methods. 

In  the  area  of  breast  tissue  modeling,  simulation  methods  have  been  developed  that 
will  allow  efficient  computation  of  scattering  from  fully  three-dimensional  models  of  breast 
tissue.  A  model  for  tissue  absorption,  important  for  realistic  simulation  of  scattering  in 
breast  tissue,  has  been  implemented.  The  frequency  dependence  of  ultrasonic  scattering 
from  soft  tissue  has  been  shown  to  be  a  critical  consideration  for  aberrating  media  such  as 
the  breast.  A  new  fc-space  method  has  been  implemented,  shown  to  provide  remarkable 
accuracy  and  efficiency,  and  extended  to  include  absorbing  boundary  conditions  as  well  as 
tissue-dependent  ultrasonic  absorption.  A  three-dimensional  implementation  of  this  method 
is  already  in  place  and  is  ready  for  computations  based  on  realistic  breast  tissue  models. 

Future  work  based  on  these  breakthroughs  will  follow  the  approved  Statement  of  Work. 
Plans  for  the  next  year  of  research  include  extension  of  time-domain  inverse  scattering  to 
include  aberration  correction,  testing  of  the  new  quantitative  imaging  algorithm  using  mea¬ 
sured  scattering  data,  construction  of  realistic  breast  tissue  models  from  cross-sectional  and 
volume  data,  and  computation  of  ultrasonic  propagation  through  these  tissue  models.  In 
the  third  year  of  the  project,  definitive  aberration-corrected  reconstructions  will  be  com¬ 
puted  both  for  measured  scattering  data  from  breast-mimicking  phantoms  and  for  simulated 
scattering  data  using  detailed  breast  models.  These  results  will  be  carefully  evaluated  for 
accuracy  and  clinical  utility,  in  collaboration  with  the  named  clinical  consultants. 

The  final  outcome  of  the  successfully  completed  project  will  be  a  novel  method  for  early 
detection,  characterization,  and  treatment  monitoring  of  breast  cancer  lesions.  Results  to 
date  indicate  that  the  final  method  will  provide  image  quality  and  diagnostic  information 
greatly  superior  to  current  2D  and  3D  ultrasonic  mammography  methods  The  finally  re¬ 
sulting  method  is  expected  to  be  competitive  with  magnetic  resonance  imaging  and  x-ray 
computed  tomography  as  a  tool  for  breast  cancer  diagnosis,  while  maintaining  inherent  ad¬ 
vantages  of  ultrasound  such  as  lower  cost,  ability  to  characterize  cystic  and  solid  lesions,  and 
safe  nonionizing  radiation. 
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Abstract 

A  quantitative  ultrasonic  imaging  method  employing  time-domain  scatter¬ 
ing  data  is  presented.  This  method  provides  tomographic  images  of  medium 
properties  such  as  the  sound  speed  contrast;  these  images  are  equivalent  to 
multiple-frequency  filtered-backpropagation  reconstructions  using  all  frequen¬ 
cies  within  the  bandwidth  of  the  incident  pulse  employed.  However,  image 
synthesis  is  performed  directly  in  the  time  domain  using  coherent  combination 
of  farfield  scattered  pressure  waveforms,  delayed  and  summed  to  numerically 
focus  on  the  unknown  medium.  The  time-domain  method  is  more  efficient 
than  multiple-frequency  diffraction  tomography  methods,  and  can,  in  some 
cases,  be  more  efficient  than  single-frequency  diffraction  tomography.  Exam¬ 
ple  reconstructions,  obtained  using  synthetic  data  for  two-dimensional  and 
three-dimensional  scattering  of  wideband  pulses,  show  that  the  time-domain 
reconstruction  method  provides  image  quality  superior  to  single-frequency 
reconstructions  for  objects  of  size  and  contrast  relevant  to  medical  imaging 
problems  such  as  ultrasonic  mammography.  The  present  method  is  closely  re¬ 
lated  to  existing  synthetic-aperture  imaging  methods  such  as  those  employed 
in  clinical  ultrasound  scanners.  Thus,  the  new  method  can  be  extended  to  in¬ 
corporate  available  image-enhancement  techniques  such  as  time-gain  compen¬ 
sation  to  correct  for  medium  absorption  and  aberration  correction  methods 
to  reduce  error  associated  with  weak  scattering  approximations. 

43.20.Fn,  43.60.Rw,  43.80.Vj,  43.20.Px 
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INTRODUCTION 

Quantitative  imaging  of  tissue  properties  is  a  potentially  useful  technique  for  diagnosis 
of  cancer  and  other  pathological  conditions.  Inverse  scattering  methods  such  as  diffraction 
tomography  can  provide  quantitative  reconstruction  of  tissue  properties  including  sound 
speed,  density,  and  absorption.  However,  although  previous  inverse  scattering  methods 
have  achieved  high  resolution  and  quantitative  accuracy,  such  methods  have  not  yet  been 
incorporated  into  commercially  successful  medical  ultrasound  imaging  systems. 

Current  inverse  scattering  methods  are  lacking  in  several  respects  with  respect  to  conven¬ 
tional  B-scan  and  synthetic  aperture  imaging  techniques.  Previous  methods  of  diffraction 
tomography,  including  methods  based  on  the  Born  approximation,  the  Rytov  approxima¬ 
tion,  and  higher-order  nonlinear  approaches,  have  usually  been  based  on  single-frequency 
scattering  (e.g.,  Refs.  1-11)  while  current  diagnostic  ultrasound  scanners  employ  wideband 
time-domain  signals.  The  use  of  wideband  information  in  image  reconstruction  is  known  to 
provide  increased  point  and  contrast  resolution,12’13  both  of  which  are  important  for  medical 
diagnosis.12’14,15 

Several  approaches  have  been  used  to  incorporate  wideband  scattering  information  into 
quantitative  ultrasonic  imaging.  One  group  of  methods  employs  time-domain  tomography 
based  on  Radon-transform  relationships  that  hold  (under  the  assumption  of  weak  scattering) 
between  scattered  acoustic  fields  and  the  reflectivity  or  scattering  strength  of  the  medium. 
Pioneering  studies  in  this  area16-18  employed  measurements  of  reflectivity  in  pulse-echo 
mode,  while  later  studies  have  incorporated  aberration  correction19,20  and  multiple-angle 
scattering  measurements.21-23  A  limitation  of  these  methods,  however,  is  that  the  Radon 
transform  relationship  strictly  holds  only  when  the  medium  is  insonified  by  an  impulsive 
(infinite  bandwidth)  wave.  When  pulses  of  finite  bandwidth  are  employed,  image  quality 
can  degrade  significantly.24 

A  number  of  linear  and  nonlinear  diffraction  tomography  methods  have  been  imple¬ 
mented  using  scattering  data  for  a  number  of  discrete  frequencies  (e.g.,  Refs.  25-32).  Al¬ 
though  use  of  multiple-frequency  data  provides  improvements  in  image  quality,  computa¬ 
tional  requirements  for  multiple-frequency  imaging  are  typically  large  because  the  computa¬ 
tional  cost  is  proportional  to  the  number  of  frequencies  employed.  To  achieve  image  quality 
competitive  with  present  diagnostic  scanners,  together  with  quantitative  imaging  of  tissue 
properties,  present  frequency-domain  methods  may  require  solution  of  the  inverse  scatter- 
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ing  problem  for  many  frequencies  within  the  bandwidth  of  the  transducer  employed.  This 
approach  thus  demands  a  high  computational  cost,  so  that  high-quality  real-time  imaging 
may  not  be  presently  feasible  using  current  frequency-domain  inverse  scattering  methods. 

Very  few  previous  workers  have  investigated  direct  use  of  time-domain  waveform  data  for 
inverse  scattering  methods  analogous  to  frequency-domain  diffraction  tomography.  Several 
methods33,34  have  used  frequency  decomposition  of  scattered  pulses  to  construct  a  wide¬ 
band  estimate  of  the  spatial  Fourier  transform  of  an  unknown  medium;  after  appropriate 
averaging  and  interpolation,  this  transform  can  be  inverted  to  obtain  a  wideband  Born 
reconstruction  of  the  medium.  A  study  reported  in  Ref.  35  has  showed  that  broadband 
synthetic  aperture  imaging  using  linear  arrays  is  closely  related  to  inverse  scattering  using 
filtered  backpropagation.  Another  method,  suggested  in  Ref.  36,  provides  a  time-domain 
reconstruction  algorithm  that  employs  filtered  backpropagation  of  scattered  fields  measured 
on  a  circular  boundary.  However,  the  time  domain  reconstruction  formula  of  Ref.  36  yields 
reconstructions  that  are  less  general  than  multiple-frequency  reconstructions  obtained  using 
the  same  signal  bandwidth. 

The  present  paper  offers  a  new  approach  to  wideband  quantitative  imaging:  a  time- 
domain  inverse  scattering  method  that  overcomes  some  of  the  limitations  of  previous 
frequency-domain  and  time-domain  quantitative  imaging  methods.  The  new  method  pro¬ 
vides  tomographic  reconstructions  of  unknown  scattering  media  using  the  entire  available 
bandwidth  of  the  signals  employed.  Reconstructions  are  performed  using  scattering  data 
measured  on  a  surface  surrounding  the  region  of  interest,  so  that  the  method  is  well  suited 
to  ultrasonic  mammography.  The  reconstruction  algorithm  is  derived  as  a  simple  delay-and- 
sum  formula  similar  to  synthetic-aperture  algorithms  employed  in  conventional  clinical  scan¬ 
ners.  However,  unlike  current  clinical  scanners,  the  present  method  can  provide  quantitative 
images  of  tissue  properties  such  as  the  spatially-dependent  sound  speed.  Reconstructions  ob¬ 
tained  in  this  manner  are  equivalent  to  reconstructions  obtained  by  combining  conventional 
frequency-domain  diffraction  tomography  reconstructions  for  all  frequencies  within  the  sig¬ 
nal  bandwidth  of  interest.  The  current  method,  however,  can  be  even  more  efficient  than 
single-frequency  diffraction  tomography.  The  method  is  applicable  both  to  two-dimensional 
and  three-dimensional  image  reconstruction.  The  direct  time-domain  nature  of  the  recon¬ 
struction  algorithm  allows  straightforward  incorporation  of  depth-  and  frequency-  dependent 
amplitude  correction  to  compensate  for  medium  absorption  as  well  as  aberration  correction 
methods  to  overcome  limits  set  by  the  Born  approximation. 
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I.  THEORY 


A.  The  time-domain  reconstruction  algorithm 


An  inverse  scattering  algorithm,  applicable  to  quantitative  imaging  of  tissue  and  other 
inhomogeneous  media,  is  derived  below.  For  simplicity  of  derivation,  the  medium  is  modeled 
as  a  fluid  medium  defined  by  the  sound  speed  variation 


(1) 


where  c0  is  a  background  sound  speed  and  c(r)  is  the  spatially-dependent  sound  speed 
defined  at  all  points  r.  For  the  scope  of  the  initial  derivation,  the  medium  is  assumed 
to  have  constant  density,  no  absorption,  and  weak  scattering  characteristics;  extensions  to 
the  reconstruction  algorithm  that  overcome  these  limiting  assumptions  are  discussed  in  the 
following  section. 

For  the  model  of  the  scattering  medium  represented  by  Eq.  (1),  the  time-domain  scattered 
acoustic  pressure  ps{r,t)  obeys  the  wave  equation37 


where  p{ r,  t)  is  the  total  acoustic  pressure  in  the  medium. 

The  scattering  configuration  considered  here  is  sketched  in  Fig.  1.  The  medium  is  sub¬ 
jected  to  a  pulsatile  plane  wave  propagating  in  the  direction  of  the  unit  vector  a, 


Pinc(r,  ot,  t)  =  f{t-  r  •  a/co), 


(3) 


where  /  is  the  time-domain  waveform  and  c0  is  the  background  sound  speed.  The  scattered 
wavefield  ps(0,  a,  t)  is  measured  at  a  fixed  radius  R  in  the  far  field,  where  6  corresponds  to 
the  direction  vector  of  a  receiving  transducer  element.  (Alternatively,  if  scattering  measure¬ 
ments  are  made  in  the  near  field,  the  farfield  acoustic  pressure  can  be  computed  using  exact 
transforms  that  represent  propagation  through  a  homogeneous  medium.25) 

A  general  time-domain  solution  for  the  wave  equation  (2),  valid  for  two-dimensional  (2D) 
or  three-dimensional  (3D)  scattering,  is  then 

/°° 

ps(0,a,Lu)e~luitdu,  (4) 

-OO 

where  ps(6 ,  a,u>)  is  a  single  frequency  component  of  the  scattered  wavefield, 
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1  r°° 

ps(0,a,u)  =  —  ps(0,a,t)eluJtdt, 

Z7T  J  —  OO 

given  exactly  by37 

ps(0,a,u))  =  k2f(co)  Jgo(R0  -  r0,u)y(r0)p{r,a,u)  dV0. 


(5) 

(6) 


In  Eq.  (6),  k  is  the  wavenumber  uj/cq  and  p( r0,  a,  u )  is  the  total  acoustic  pressure  associated 
with  the  unit-amplitude  incident  plane  wave  elfcar°.  The  integral  in  Eq.  (6)  is  taken  over 
the  entire  support  of  7  in  5ft2  for  2D  scattering  or  in  5ft3  for  3D  scattering.  The  free-space 
Green’s  function,  represented  by  Go  in  Eq.  (6),  is38 


Go(r,u;)  =  ^  H^ikr)  for  2D  scattering  and 

eikr 

G0( t,uj)  =  - —  for  3D  scattering,  (7) 

47rr 

where  is  the  zeroth-order  Hankel  function  of  the  first  kind  and  r  is  the  magnitude  of 
the  vector  r. 

The  farfield  scattered  pressure,  when  specified  for  all  incident- wave  directions  a,  mea¬ 
surement  directions  0,  and  times  t ,  comprises  the  data  set  to  be  used  for  reconstruction  of 
the  unknown  medium  7(r).  The  inverse  scattering  problem  is  to  reconstruct  7(1*)  using  this 
measured  data.  The  starting  point  for  the  present  time-domain  inverse  scattering  method 
is  conventional  single-frequency  diffraction  tomography,  which  can  be  performed  using  any 
frequency  component  of  the  measured  farfield  scattering. 

Under  the  assumption  of  weak  scattering,  one  can  make  the  Born  approximation,  in 
which  the  total  pressure  p(cc,  u>)  in  Eq.  (6)  is  replaced  by  the  plane  wave  elkr  a.  For  scattering 
measurements  made  at  a  radius  R  in  the  far  field,  linearized  inversion  of  Eq.  (6)  for  any 
frequency  component  can  then  be  performed  using  filtered  backpropagation, 4,25,39,40  i.e., 

7a(r,w)  =  -  //  m°‘)pA9,a,“)e‘k(e-a)" dSedS,,  (8) 

/M  JJ 

where 


$(0,a) 

HO,  a) 


|  sin(#  —  ck)|  in  2D,  and 
\0  —  a |  in  3D. 


(9) 


Each  surface  integral  in  Eq.  (8)  is  performed  over  the  entire  measurement  circle  for  the  2D 
case  and  over  the  entire  measurement  sphere  for  the  3D  case. 
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Equation  (8)  provides  an  exact  solution  to  the  linearized  inverse  scattering  problem 
for  a  single  frequency  component  of  the  scattered  wavefield  ps(0,et,t).  Instead  of  seeking  a 
solution  for  7(1*)  that  exactly  satisfies  the  linearized  time-domain  inverse  scattering  problem, 
the  method  outlined  here  constructs  a  “compounded”  7(r)  by  integrating  single-frequency 
reconstructions  7(1*,  u>)  over  a  range  of  frequencies  c 0.  A  generalized  formula  for  this  approach 


can  be  written 


7m  (r) 


/o0°g(^)7B(r,a;)da; 
/o°°  9(a>)  du, 


where  g(u)  is  an  appropriate  frequency-dependent  weighting  function. 
Using  Eq.  (8),  and  making  the  definition 


roo 

N  =  2  g(u>)  du>, 

Jo 


Eq.  (10)  can  be  written  in  the  form 


7mM  =  |  l°°m  *(».<*)  0.(0.  c,  dS^dSedn,.  (12) 

If  the  frequency  weight  g(u>)  is  now  specified  to  incorporate  the  incident-pulse  spectrum 
f(oo)  and  to  compensate  for  the  frequency-  and  dimension-dependent  coefficient 

g(u)  =  (13) 


Eq.  (12)  reduces  to  the  form 

7A/W  =  ^IJ  <He,o)  jT  ps(9,a,w)e-“lR+<«-8'']rfwdSo(iS».  (14) 

Equation  (14)  can  be  rewritten  using  the  definition  oips(0,  a,u>)  from  Eq.  5  to  yield  the 
time-domain  expression 

7mM  =  ^  //«■(«,  a)  L  p,  (»,  a,  R/c„  +  (°  '  r)  dS„  dS ,,  (15) 

where  L  denotes  the  linear  operator 

roo  A 

L  hp(t)]  =  2  /  i[)(u))e^'lut  dee  (16) 

Jo 

and  ip(co)  is  the  Fourier  transform  of  ip(t)  using  the  definition  from  Eq.  (5). 

Using  the  conjugate  symmetry  of  ^(oo)  [i i.e .,  ip(0,  a, co)  —  ip*(9,  a,  —u>)  for  any  real 
the  real  part  of  L[ip(t)]  is  shown  to  be  simply  Similarly,  using  the  convolution  theorem 
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> 


as  well  as  the  conjugate  symmetry  of  ip(t),  the  imaginary  part  of  L[ip(t)]  is  seen  to  be  an 
inverse  Hilbert  transform41  of 


Im 


L  [V-(*)l 


fA  dr  =  H-1  [«()]. 


(17) 


This  transform,  also  known  as  a  quadrature  filter,  applies  a  phase  shift  of  7t/2  to  each 
frequency  component  of  the  input  signal. 

Thus,  the  time-domain  reconstruction  formula  can  finally  be  written 

7m (r)  =  jj  JJ  $(0,  Qf)  ( Ps{0 ,  a,  r)  +  iH-1  \ps(0,  a,  r)])  dSa  dSe ,  (18) 


where 

r  =  fl/c  +  AjAtl.  (19) 

Co 

The  direction-dependent  weight  $(0,  a),  which  is  the  same  as  the  “filter”  employed  in 
single-frequency  filtered  backpropagation,  is  given  for  the  2D  and  3D  cases  by  Eq.  (9). 

Equation  18  is  notable  in  several  respects.  First,  it  provides  a  linearized  reconstruction 
that  employs  scattering  information  from  the  entire  signal  bandwidth  without  any  frequency 
decomposition  of  the  scattered  wavefield.  Second,  the  delay  term  r  corresponds  exactly 
to  the  delay  required  to  construct  a  focus  at  the  point  r  by  delaying  and  summing  the 
scattered  wavefield  ps(6,  ct,  t)  for  all  measurement  directions  6  and  incident-wave  directions 
a.  Thus,  the  time-domain  reconstruction  formula  given  by  Eq.  (18)  can  be  regarded  as  a 
quantitative  generalization  of  confocal  time-domain  synthetic  aperture  imaging,  in  which 
signals  are  synthetically  delayed  and  summed  for  each  transmit/receive  pair  to  focus  at  the 
image  point  of  interest.35,42-45 

A  reconstruction  formula  similar  to,  although  less  general  than,  Eq.  18  was  independently 
derived  in  Ref.  36  for  the  two-dimensional  inverse  scattering  problem.  In  view  of  the  present 
derivation,  the  method  of  “probing  by  plane  pulses”  in  Ref.  36  can  be  regarded  to  yield 
a  frequency-compounded  reconstruction  of  Re  [7m (r)],  while  the  present  method  yields  the 
complex  function  7m(f).  In  Ref.  36,  this  method  was  proposed  as  a  more  convenient  way 
to  implement  narrow-band  diffraction  tomography.  However,  the  numerical  results  given 
below  show  that  the  reconstruction  formula  of  Eq.  18,  when  directly  implemented  using 
wideband  signals,  provides  considerable  improvement  in  image  quality  over  narrow- band 
reconstructions. 
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Reconstructions  using  Eq.  (18)  can  be  performed  using  any  pulse  waveform.  However, 
the  frequency  compounding  defined  by  Eq.  (10)  is  most  straightforwardly  interpreted  if  the 
frequency  weight  g(u)  has  a  phase  that  is  independent  of  frequency.  This  criterion  can  be 
met,  for  instance,  if  the  incident  pulse  waveform  f(t)  is  even  in  time, 

m  =  /(-«),  (20) 

so  that  f(u)  is  purely  real.  (Similarly,  if  the  incident  pulse  waveform  is  odd  in  time,  f(u)  is 
purely  imaginary  and  Eq.  (18)  can  still  be  employed.) 

However,  supposition  of  a  frequency-independent  phase  for  /( co)  does  not  result  in  any 
loss  of  generality.  For  any  linear-phase  signal,  such  that  the  Fourier  transform  has  the  form 

/»  =  |/H|e"<,  w>0,  (21) 

an  additional  delay  term  of  magnitude  £  can  be  applied  to  all  scattered  signals  to  obtain 
the  signals  associated  with  the  purely-real  spectrum  |/(w)|.  In  general,  the  scattered  field 
associated  with  a  desired  waveform  f(t)  can  be  determined  for  an  arbitrary  waveform  u(t) 
from  the  deconvolution  operation 

[p,(@,  a,  «)],<„  =F"‘  [?»(#,«, *)]„,,,.  (22) 

For  stable  deconvolution  using  Eq.  (22),  the  desired  f(u)  should  not  have  significant  fre¬ 
quency  components  outside  the  bandwidth  of  u(u). 

B.  Extensions  to  the  reconstruction  algorithm 

For  large  tissue  structures  at  high  ultrasonic  frequencies,  weak  scattering  approximations 
such  as  the  Born  approximation  are  of  limited  validity.  Thus,  for  problems  of  interest  to 
medical  ultrasound  imaging,  reconstructed  image  quality  can  be  improved  by  aberration 
correction  methods  that  incorporate  higher-order  scattering  and  propagation  effects.  The 
present  time-domain  reconstruction  formula  (18)  provides  a  natural  framework  for  quanti¬ 
tative  imaging  with  aberration  correction.  In  general,  if  the  background  medium  is  known 
or  can  be  estimated,  the  received  scattered  signals  can  be  processed  to  provide  an  estimate 
of  the  scattered  field  that  would  be  measured  for  the  same  scatterer  within  a  homogeneous 
background  medium.  This  approach  essentially  removes  higher-order  scattering  effects  from 
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the  measured  far  field  scattering,  so  that  a  Born  inversion  can  be  performed  on  the  modified 
data;  similar  processes  occur  implicitly  in  many  nonlinear  inverse  scattering  methods.46 

For  example,  a  simple  implementation  of  aberration  correction  can  be  derived  if  one 
makes  the  assumption  that  background  inhomogeneities  result  only  in  cumulative  delays 
(or  advances)  of  the  incident  and  scattered  wavefronts.  This  crude  model  does  not  include 
many  propagation  and  scattering  effects  important  to  ultrasonic  aberration,  but  has  been 
shown  to  provide  a  reasonable  first  approximation  of  local  delays  in  wavefronts  propagating 
through  large-scale  tissue  models.47,48  Given  this  approximation,  the  total  delay  for  an  angle 
<fi  and  a  point  position  r  is  given  by 

6t(<P,  r)  =  f  c(£)_1  df  -  —  (23) 

H  Co 

where  the  integral  is  performed  along  the  line  that  joins  the  spatial  points  r  and  Re/), 
Aberration-corrected  reconstructions  can  then  be  performed  using  Eq.  (18)  with  r  replaced 
by  the  corrected  delay  term 

r ->• -R/c0  +  — — —  -  +  5r(a,r)  +  Sr(d,r).  (24) 

Co 

Improved  approximations  could  be  obtained  by  application  of  the  delay  function  Sr(4>,  r) 
after  numerical  backpropagation  of  the  farfield  scattered  wavefronts  through  a  homogeneous 
medium49,50  or  by  compensation  for  both  delay  and  amplitude  variations.51’52  More  general, 
although  much  more  computationally  expensive,  aberration  correction  could  also  be  per¬ 
formed  by  synthetic  focusing  using  full-wave  numerical  computation  of  acoustic  fields  within 
an  estimated  realization  of  the  unknown  medium.  This  method  has  been  implemented, 
within  the  context  of  a  frequency-domain  diffraction  tomography  method,  in  Ref.  30. 

The  present  imaging  method  has  been  derived  using  simplifying  assumptions  including 
zero  absorption  and  constant  density  for  the  scattering  medium.  However,  these  assumptions 
do  not  substantially  restrict  the  validity  of  the  method.  For  example,  the  effect  of  absorp¬ 
tion  can  be  reduced  using  time-gain  compensation,  with  or  without  frequency-dependent 
corrections,53  of  received  scattered  signals  for  each  transmit/receive  pair.  Such  time-gain 
compensation  could  be  performed  either  using  an  estimated  bulk  attenuation  for  the  medium 
(as  with  current  clinical  ultrasound  scanners),  or  by  implementation  of  an  adaptive  attenu¬ 
ation  model  in  a  manner  similar  to  the  time-shift  compensation  scheme  discussed  above. 

Inclusion  of  density  variations  as  well  as  sound  speed  variations  adds  additional  complica¬ 
tion  to  the  time-domain  diffraction  tomography  algorithm  derived  here.  For  single-frequency 
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diffraction  tomography  in  the  presence  of  sound  speed  and  density  variations,  the  quantity 
7b  (r,  u>)  reconstructed  by  Eq.  8  can  be  shown54  to  provide  an  estimate  of  a  physical  quantity 
that  depends  both  on  sound  speed  variations  and  density  variations.  In  the  notation  used 
here,  this  quantity  can  be  written 

7'(r)  =  7(r)  -  7(r)7p(r)  +  ^V27P(r),  (25) 

where  the  density  variation  is  defined  7p  =  1— p0/p(r).  Thus,  for  time-domain  reconstructions 
of  media  with  density  variations,  the  reconstruction  formula  of  Eq.  18  will  provide  the 
estimate 

7 w(r)  ~  7(r)  -  7(r)7p(r)  +  V27p(r),  (26) 

where  ko  is  the  wavenumber  corresponding  to  the  center  frequency  of  the  pulse  employed.  For 
media  such  as  human  tissue,  where  density  variations  are  fairly  small  and  abrupt  density 
transitions  are  rare,  the  last  two  terms  of  Eq.  (26)  are  small  compared  to  7(r),  so  that 
the  reconstruction  algorithm  derived  above  can  still  be  regarded  to  provide  an  image  of 
the  sound-speed  variation  function  7(r).  However,  if  desired,  a  reconstruction  employing 
pulses  with  two  distinct  center  frequencies  could  allow  separation  of  sound  speed  and  density 
variations  by  techniques  similar  to  those  described  in  Refs.  25  or  54. 
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II.  COMPUTATIONAL  METHODS 

The  time-domain  inverse  scattering  method  described  above  has  been  tested  with  2D  and 
3D  synthetic  data  prepared  using  three  numerical  methods:  a  Born  approximation  method 
for  point  scatterers  and  3D  slabs,  an  exact  series  solution  for  cylindrical  inhomogeneities, 
and  a  A;-space  method  for  arbitrary  2D  inhomogeneous  media. 

The  time-domain  waveform  employed  for  all  the  computations  reported  here  was 


f(t)  =  cos(u>0t)  e 


-t2/(2(72) 


where  u>0  =  27t/0  for  a  center  frequency  of  /o  and  o  is  the  temporal  Gaussian  parameter. 
This  waveform  has  the  real,  even  Fourier  transform 


/M  = 


3—cr2  (w—it>o)2/2  .  — ct2(w+u>o)2/2 


Values  used  for  the  compurations  reported  here  were  /o  =  2.5  MHz  and  o  =  0.25  /xs,  so  that 
the  —6  dB  bandwidth  of  the  signal  was  1.5  MHz.  These  parameters  correspond  closely  to 
those  of  an  existing  2048-element  ring  transducer.55 

For  the  case  of  point  scatterers,  the  contrast  function  7  was  assumed  to  take  the  form 

M 

7(r)  ~Ti)-  (29) 

1 

Using  the  farfield  form  of  the  2D  Green’s  function  and  neglecting  multiple  scattering,  Eq.  (6) 
for  the  scattered  far  field  can  be  rewritten  as 


ps(0,OL,Uj) 


8irkR 


.  Jk(a-0)'rj 


for  each  frequency  component  of  interest.  Time-domain  waveforms  were  synthesized  by  using 
Eq.  (30)  for  each  frequency  with  f(u>)  >  10-3  and  inverting  the  frequency-domain  scattered 
wavefield  by  a  fast  Fourier  transform  (FFT)  implementation  of  Eq.  (4).  The  temporal 
sampling  rate  employed  was  10  MHz.  An  analogous  formula,  with  a  different  multiplicative 
constant,  was  also  employed  for  the  3D  case. 

The  Born  approximation  was  also  used  to  compute  three-dimensional  scattering  for  slab¬ 
shaped  objects  defined  by  the  equation 

q(r)  =  70 H(ax  -  \x\)H(ay  -  \y\)H(az  -  |s|).  (31) 
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x 


(32) 


For  this  object,  the  linearized  forward  problem  can  be  solved  analytically.  Under  the  Born 
approximation,  the  frequency-domain  scattered  far  field  has  the  form 

Ps(@,  OC,  iS)  =  2  /(w)  '70  CLXQ’yQ'z  c  /(7tJ2) 

sin [kLx  (a-0)  •  ex]  sin [kLy  (a-0)  ■  ej  sin[fcLz  (a-0)  •  ez] 
kLx(a-0)-ex  kLy\a-0)-ey  kLz{a-0)-ez  ’ 
where  ex,  ey,  and  e2  represent  unit  vectors  in  the  x,  y,  and  z  directions.  The  time  domain 
scattered  pressure  ps(0,a,t)  is  obtained,  as  for  the  point  scatterer  case  described  above, 
by  inverse  transformation  of  the  frequency-domain  wavefield  for  all  frequencies  within  the 
bandwidth  of  interest. 

For  2D  cylindrical  inhomogeneities,  an  analogous  procedure  was  followed,  except  that 
the  frequency-domain  scattered  wavefield  ps(0,a,u)  was  computed  using  an  exact  series 
solution37  for  each  frequency  component  of  interest.  In  implementation  of  the  series  solution, 
summations  were  truncated  when  the  magnitude  of  a  single  coefficient  dropped  below  10~12 
times  the  sum  of  all  coefficients. 

Solutions  were  also  obtained  for  arbitrary  2D  inhomogenous  media  using  a  time-domain 
Axspace  method.56  Grid  sizes  of  256  x  256  points,  a  spatial  step  of  0.0833  mm,  and  a  time 
step  of  0.02734  ps  were  employed.  Scattered  acoustic  pressure  signals  on  a  circle  of  virtual 
receivers  were  recorded  at  a  sampling  rate  of  9.144  MHz.  The  receiver  circle,  which  had 
a  radius  of  3.0  mm  in  these  computations,  completely  contained  the  inhomogeneities  used. 
Farfield  waveforms  were  computed  by  Fourier  transforming  the  time-domain  waveforms  on 
the  nearfield  measurement  circle,  transforming  these  to  farfield  waveforms  for  each  frequency 
using  a  numerically  exact  transformation  method,25  and  performing  inverse  Fourier  transfor¬ 
mation  to  yield  time-domain  farfield  waveforms.  All  forward  and  inverse  temporal  Fourier 
transforms,  as  well  as  angular  transforms  occurring  in  the  nearfield-farfield  transformation,25 
were  performed  by  FFT. 

The  time-domain  imaging  method  was  directly  implemented  using  Eq.  (18),  evaluated  us¬ 
ing  straightforward  numerical  integration  over  all  transmit  and  receive  directions  employed. 
The  reconstruction  formula  employed  can  be  explicitly  written  as 

1m{t)  =  /  /  |  sin(o!  -  0)1  (ps(0,a,r)  +  iH-1  [ps(0,  a,r)])  dadd, 

i\2d  Jo  Jo  v 


N2D 

R/ Co  + 


(cos  a  —  cos  9)  -  x  +  (sin  a  —  sin  9)  -y 


c0 


(33) 


for  the  2D  case,  where  a  and  6  are  the  angles  corresponding  to  the  direction  vectors  a  and 
6,  and  as 
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\  p2tT  P7T  r27T  P7T  .  . 

1m{v)  =  WdL  Jo  Jo  Jo  |a  "  0|  {Ps^^r)  +  \ps(0,a,r)}) 


T  =  R/C0  + 


x  sin($a)  sin($@)  dOa  d$e  dQo  , 
(a  —  0)  ■  r 


co 


a  —  6  =  (cos  0a  sin  <Fa  —  cos  0#  sin  <3>0)  •  ex  +  (sin  ©a  sin  —  sin  0#  sin  $0)  ■  ey 

+  (cos  -  cos  $0)  •  ez  (34) 


for  the  3D  case,  where  0a  and  $a  are  direction  angles  for  the  incident-wave  direction  a 
and  00  and  $0  are  direction  angles  for  the  measurement  direction  6.  For  each  case,  the 
normalization  factor  N  was  determined  from  Eq.  11  with  g(uj)  =  / (u) / and  h(uj) 
given  by  Eq.  9.  Before  evaluation  of  the  argument  r  for  each  signal,  the  time-domain 
waveforms  were  resampled  at  a  sampling  rate  of  16  times  the  original  rate.  This  resampling 
was  performed  using  FFT-based  Fourier  interpolation.  The  inverse  Hilbert  transform  was 
performed  for  each  signal  using  an  FFT  implementation  of  Eq.  (16).  Values  of  the  pressure 
signals  at  the  time  r  were  then  determined  using  linear  interpolation  between  samples  of  the 
resampled  waveforms.  The  integrals  of  Eqs.  (33)  and  (34)  were  implemented  using  discrete 
summation  over  all  transmission  and  measurement  directions  employed. 

Computations  were  also  performed  using  the  time-domain  diffraction  tomography  algo¬ 
rithm  for  limited-aperture  data.  For  these  reconstructions,  the  integrals  of  Eq.  (33)  were 
evaluated  only  for  angles  corresponding  to  transmitters  and  receivers  within  a  specified 
aperture  of  angular  width  </>ap,  i.e., 


|c>j|  ^  ^ap/2, 
\0  -  A  <  ^ap/2, 


(35) 


Use  of  a  small  value  for  ^ap  corresponds  to  use  of  a  small  aperture  in  pulse-echo  mode. 
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III.  NUMERICAL  RESULTS 

Two-dimensional  and  three-dimensional  point-spread  functions  (PSF)  for  the  present 
time-domain  diffraction  tomography  method  are  illustrated  in  Fig.  2.  The  time-domain  re¬ 
constructions  shown  here,  like  the  other  time-domain  reconstructions  shown  in  this  paper, 
were  obtained  using  a  incident  pulse  of  center  frequency  2.5  MHz  and  a  Gaussian  envelope 
corresponding  to  a  -6  dB  bandwidth  of  1.5  MHz.  Point-spread  functions  were  determined 
by  reconstructing  a  point  scatterer  located  at  the  origin.  For  the  2D  case,  in  which  the  point 
scatterer  can  be  regarded  as  a  thin  wire,  synthetic  scattering  data  was  obtained  using  the 
Born  approximation  method  outlined  above  for  16  incident-wave  directions  and  64  measure¬ 
ment  directions.  The  3D  time-domain  reconstruction  was  obtained  using  Born  data  for  72 
incident-wave  directions  and  288  measurement  directions,  each  evenly  spaced  on  a  rectan¬ 
gular  grid  defined  by  the  angles  0  and  $.  For  comparison,  analogous  point-spread  functions 
are  also  shown  for  standard  frequency-domain  diffraction  tomography  reconstructions  with 
single-frequency  (2.5  MHz)  data. 

For  the  2D  case  illustrated  in  Fig.  2,  the  time-domain  reconstruction  has  a  slightly 
narrower  peak,  indicating  that  point  resolution  has  been  slightly  improved  by  the  increased 
bandwidth  employed  in  the  time  domain  method.  More  significantly,  sidelobes  of  the  time- 
domain  PSF  are  dramatically  smaller  than  those  for  the  frequency-domain  PSF,  so  that 
contrast  resolution  for  time-domain  diffraction  tomography  is  seen  to  be  much  higher  than 
for  single-frequency  diffraction  tomography.  For  the  3D  case,  the  time-domain  reconstruction 
shows  a  much  more  dramatic  improvement  over  the  frequency-domain  reconstruction.  In  this 
case,  both  the  point  resolution  and  the  contrast  resolution  are  significantly  higher  for  time- 
domain  reconstruction  than  for  single-frequency  reconstruction.  Furthermore,  a  comparison 
of  the  PSF’s  for  2D  and  3D  time-domain  reconstruction  indicates  that  much  higher  image 
quality  is  achievable  for  3D  time-domain  imaging  than  for  the  2D  case.  This  increase  in 
image  quality  suggests  that  the  time-domain  diffraction  tomography  method  proposed  here 
may  benefit  from  the  overdetermined  nature  of  the  general  wideband  3D  inverse  scattering 
problem.57,58 

The  effect  of  transmit  and  receive  aperture  characteristics  on  image  quality  is  illustrated 
in  Fig.  3,  which  shows  the  point-spread  function  for  a  number  of  aperture  configurations. 
All  of  these  reconstructions  employed  64  measurement  directions,  which  were  found  to  ad¬ 
equately  sample  the  scattered  far  field  for  a  single  transmit  direction.  Figure  3(a)  shows 
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the  point-spread  function  for  reconstructions  obtained  using  1,  4,  8,  and  16  incident-wave 
directions.  The  point  scatterer  is  clearly  imaged  even  for  the  reconstruction  using  one 
incident- wave  direction.  Optimal  image  quality  (indistinguishable  from  reconstructions  with 
64  incident- wave  directions)  is  obtained  for  16  incident- wave  directions,  so  that  scattering 
data  obtained  using  one  incident-wave  direction  for  each  group  of  four  measurement  direc¬ 
tions  appears  to  be  sufficient  for  the  present  reconstruction  method.  The  effect  of  limited 
view  range  on  the  point  spread  function  is  illustrated  in  Fig.  3.  Limitation  of  the  transmit 
and  recieve  apertures  to  angles  near  the  backscatter  direction  (aperture  size  n/2)  results 
in  an  image  that  resembles  a  conventional  B-scan.  Use  of  an  aperture  corresponding  to 
pulse-echo  mode  in  the  large-aperture  limit  (aperture  size  7 r)  yields  higher  resolution  in  all 
directions.  Using  three  fourths  of  a  circular  aperture  (size  3n/2)  yields  image  quality  close 
to  that  for  the  full  aperture  (2n)  case. 

Reconstructions  performed  using  exact  solutions  for  scattering  from  cylindrical  inho¬ 
mogeneities  provide  a  straightforward  means  to  assess  the  accuracy  of  the  time-domain 
scattering  method  for  a  range  of  object  sizes  and  contrasts.  Cross  sections  of  time-domain 
and  single-frequency  reconstructions,  plotted  in  Fig.  4,  show  the  relative  accuracy  of  each 
reconstruction  method  for  a  cylinder  of  1  mm  radius  and  purely  real  contrast  ranging  from 
7  =  0.02  to  7  =  0.08.  For  the  synthetic  scattering  data  in  each  case,  96  measurement  di¬ 
rections  and  24  incident-wave  directions  were  employed.  The  time-domain  reconstructions 
show  improvement  over  the  single-frequency  reconstructions  both  in  improved  contrast  reso¬ 
lution  (smaller  sidelobes  outside  the  support  of  the  cylinder)  and  in  decreased  ripple  (Gibbs 
phenomenon)  artifacts  within  the  support  of  the  cylinder.  However,  for  increasing  con¬ 
trast  values,  both  methods  show  similar  increases  in  phase  error,  as  indicated  by  increased 
imaginary  parts  of  the  reconstructed  contrast. 

Further  testing  of  accuracy  for  the  time-domain  reconstruction  method  is  shown  in  Fig.  5. 
Real  parts  of  time-domain  reconstructions  are  shown  for  cylinder  radii  between  1  mm  and 
4  mm  with  contrasts  between  7  =  0.02  and  7  =  0.08.  The  number  of  measurement  directions 
employed  for  the  synthetic  scattering  data  was  96  for  the  1  mm  radius  cylinders,  192  for  the 
2  mm  cylinders,  288  for  the  3  mm  cylinders,  and  384  for  the  4  mm  cylinders.  In  each  case, 
four  incident-wave  directions  per  measurement  direction  were  used.  Using  the  wavenumber 
k0  =  10.472  rad/mm  corresponding  to  the  center  frequency  of  2.5  MHz  and  a  sound  speed  of 
1.5  mm/jus,  the  reconstructions  shown  in  Fig.  5  indicate  that  the  interior  of  the  cylinder  is 
imaged  satisfactorily  for  the  approximate  range  kax  7  <  2.5.  This  result  is  consistent  with 
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a  previous  study  of  single-frequency  diffraction  tomography,  which  showed  that  adequate 
Born  reconstructions  of  cylinders  were  obtained  for  the  parameter  range  kax 7  <  2. 2. 59 

Reconstructions  for  several  scattering  objects  without  special  symmetry  are  shown  in 
Fig.  6.  All  of  these  reconstructions  were  performed  using  synthetic  data  produced  by  the  k- 
space  method  described  in  Ref.  56.  Synthetic  scattering  data  were  computed  for  64  incident- 
wave  directions  and  256  measurement  directions  in  each  case.  The  first  panel  shows  a 
reconstruction  of  a  cylinder  of  radius  2.5  mm  and  contrast  7  =  —0.0295  with  an  internal 
cylinder  of  radius  0.2  mm  and  contrast  7  =  0.0632.  The  second  panel  shows  a  reconstruction 
of  a  2.5  mm- radius  cylinder  with  random  internal  structure.  The  third  reconstruction  shown 
employed  a  portion  of  a  chest  wall  tissue  map  from  Ref.  60.  In  this  case,  the  synthetic 
data  was  obtained  using  a  tissue  model60  that  incorporates  both  sound  speed  and  density 
variations,  so  that  the  reconstructed  quantity  is  given  by  Eq.  (26). 

The  real  part  of  each  reconstruction  shows  good  image  quality,  with  high  resolution  and 
very  little  evidence  of  artifacts.  Particularly  notable  is  the  accurately  detailed  imaging  of 
internal  structure  for  the  random  cylinder  and  the  chest  wall  cross  section.  As  expected, 
the  density  variations  present  in  the  chest  wall  cross  section  have  not  greatly  affected  the 
image  appearance;  there  is,  however,  a  slight  edge  enhancement,  associated  with  the  Lapla- 
cian  term  in  Eq.  (26),  at  boundaries  between  tissue  regions.  Also  notable  is  the  nearly- 
complete  absence  of  any  artifacts  outside  the  scatterer  in  each  case;  this  result  indicates 
that  high  contrast  resolution  has  been  achieved.  However,  in  each  case,  the  imaginary  part 
of  the  reconstruction  is  nonzero,  indicating  that  the  Born  approximation  is  not  fully  appli¬ 
cable.  Aberration  correction  methods  [of  which  a  simple  example  is  given  by  Eq.  (24)]  could 
substantially  reduce  this  phase  error,  as  for  multiple-frequency  diffraction  tomography  in 
Ref.  30. 

Three-dimensional  reconstructions  of  a  homogeneous  slab  are  shown  in  Fig.  7.  The 
scatterer  is  characterized  by  Eq.  (31)  with  70  =  0.01,  ax  =  0.5  mm,  ay  =  1.0  mm,  and 
az  —  1.5  mm.  Synthetic  data  was  computed  using  Eq.  34  for  288  transmit  directions  and 
1152  measurement  directions,  each  evenly  spaced  in  the  angles  $  and  0.  Signal  parameters 
were  as  for  the  examples  above,  except  that  the  initial  sampling  rate  for  the  time-domain 
signals  was  9.0  MHz.  Isosurface  renderings  of  the  real  part  of  7 m  are  shown  for  the  sur¬ 
faces  7m  =  0.0025.  Since  the  scattering  data  were  obtained  using  a  Born  approximation 
for  the  3D  case,  the  imaginary  part  of  each  reconstruction  is  identically  zero  for  both  re¬ 
constructions.  Consistent  with  the  point-spread  functions  shown  in  Fig.  2,  the  time-domain 
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reconstruction  is  much  more  accurate  than  the  single-frequency  reconstruction.  While  the 
single-frequency  reconstruction  shows  an  erroneously  rippled  surface,  the  time-domain  re¬ 
construction  is  smooth.  The  time-domain  reconstruction  is  nearly  identical  to  the  original 
object  except  for  some  rounding  of  edges  due  to  the  limited  high-frequency  content  of  the 
signal  employed. 

Since  three-dimensional  inverse  scattering  is  a  computationally  demanding  problem,  com¬ 
parison  of  computational  efficiency  for  single-frequency  and  time-domain  methods  is  of  inter¬ 
est.  For  both  reconstructions  shown  in  Fig.  7,  identical  discretizations  of  the  reconstructed 
medium  were  employed.  Both  computations  included  solution  of  the  applicable  linearized 
forward  problem  as  well  as  the  inverse  problem.  Nonetheless,  the  time-domain  method 
was  more  efficient  than  the  single-frequency  method;  the  total  CPU  time  required  on  a 
200  MHz  AMD  K6  processor  was  133.3  CPU  min  for  the  time-domain  method  and  287.4 
CPU  minutes  for  the  single-frequency  method.  This  gain  in  efficiency  was  possible  because 
the  greatest  computational  expense  occurred  in  the  “backpropagation”  of  the  signals  for 
each  reconstruction  point.  For  the  single-frequency  method,  this  step  required  evaluation  of 
complex  exponentials  for  each  transmit  direction,  measurement  direction,  and  spatial  point. 
For  the  time-domain  method,  however,  the  computationally  intensive  steps  (including  the 
forward  problem  solution  and  Fourier  interpolation  of  the  scattered  signals)  needed  only  to 
be  performed  once  for  each  transmit/receive  pair.  For  the  backpropagation  step,  performed 
at  each  point  in  the  3D  spatial  grid,  the  time-domain  reconstruction  method  required  only 
linear  interpolation  of  the  oversampled  farfield  pressure  waveforms. 
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IV.  CONCLUSIONS 

A  new  method  for  time-domain  ultrasound  diffraction  tomography  has  been  presented. 
The  method  provides  quantitative  images  of  sound  speed  variations  in  unknown  media; 
when  two  pulse  center  frequencies  are  employed,  the  method  is  also  capable  of  imaging 
density  variations.  Reconstructions  performed  using  this  method  are  equivalent  to  multiple- 
frequency  reconstructions  using  filtered  backpropagation,  but  can  be  obtained  with  much 
greater  efficiency. 

The  time-domain  reconstruction  algorithm  has  been  derived  as  a  simple  filtered  delay- 
and-sum  operation  applied  to  farfield  scattered  signals.  This  algorithm  is  closely  related  to 
time-domain  confocal  synthetic  aperture  imaging,  so  that  it  can  be  considered  a  generaliza¬ 
tion  of  imaging  algorithms  employed  in  current  clinical  instruments.  The  simplicity  of  the 
imaging  algorithm  allows  straightforward  addition  of  features  such  as  time-gain  compensa¬ 
tion  and  aberration  correction. 

Numerical  results  obtained  using  synthetic  data  for  2D  and  3D  scattering  objects  show 
that  the  time-domain  method  can  yield  significantly  higher  image  quality  (and,  in  some 
cases,  also  greater  efficiency)  than  single-frequency  diffraction  tomography.  Quantitative  re¬ 
constructions,  obtained  using  signal  parameters  comparable  to  those  for  present-day  clinical 
instruments,  show  accurate  imaging  of  objects  with  simple  deterministic  structure,  random 
internal  structure,  and  structure  based  on  a  cross-sectional  tissue  model.  The  method  is 
hoped  to  be  useful  for  diagnostic  imaging  problems  such  as  the  detection  and  characteriza¬ 
tion  of  lesions  in  ultrasonic  mammography. 
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FIGURES 


FIG.  1.  Scattering  configuration.  An  incident  pressure  pulse  of  the  form  f(t  —  a  ■  r/c)  is 
scattered  by  an  inhomogeneous  medium  and  the  time-domain  scattered  field ps{Q,  a,  t)  is  measured 
at  a  radius  R  in  the  far  field. 
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FIG.  2.  Point-spread  functions  for  time-domain  and  single-frequency  diffraction  tomography 
methods.  In  each  panel,  the  vertical  scale  corresponds  to  the  relative  amplitude  of  the  recon¬ 
structed  contrast  7(1-),  while  the  horizontal  scale  corresponds  to  number  of  wavelengths  at  the 
center  frequency,  (a)  Two-dimensional  case,  (b)  Three-dimensional  case. 
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FIG.  3.  Effect  of  aperture  characteristics  on  image  quality.  The  point-spread  function  (PSF), 
determined  as  the  real  part  of  the  reconstructed  contrast  7M(r),  is  shown  for  the  same  waveform 
parameters  as  in  Fig.  2.  Each  panel  shows  an  area  of  0.6  mm  x  0.6  mm,  corresponding  to  one  square 
wavelength  at  the  center  frequency.  In  each  case,  64  measurement  directions  were  employed, 
(a)  PSF  for  1,  4,  8,  and  16  incident- wave  directions,  (b)  PSF  for  aperture  sizes  of  7r/2,  7 r,  3n/2, 


and  2-7T  radians. 
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FIG.  4.  Cross  sections  of  reconstructed  contrast  functions  7(1-)  for  a  cylinder  of  radius  1  mm, 
using  time-domain  (TD)  and  single-frequency  (SF)  diffraction  tomography.  Waveform  parameters 
are  as  in  Fig.  1.  (a)  7  =  0.02.  (b)  7  =  0.04.  (c)  7  =  0.06.  (d)  7  =  0.08. 
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FIG.  5.  Time-domain  reconstructions  for  cylinders  of  varying  size  and  contrast.  Each  panel 
shows  the  real  part  of  the  reconstructed  contrast  Jm(v)  for  a  pulse  of  center  frequency  2.5  MHz 
and  -6  dB  bandwidth  1.5  MHz.  In  each  case,  the  ratio  of  measurement  directions  to  incident- wave 
directions  is  four.  All  images  axe  shown  on  a  linear,  bipolar  gray  scale  where  white  represents  the 
maximum  amplitude  of  |7m(f)|  and  black  represents  —1  times  the  maximum  amplitude. 
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FIG.  6.  Time-domain  reconstructions  from  full-wave  synthetic  data  for  three  arbitrary  scatter¬ 
ing  objects.  Each  row  shows  the  actual  (purely  real)  contrast  function  7  together  with  the  real 
and  imaginary  parts  of  the  reconstructed  contrast  function  7 m,  using  the  same  linear  bipolar  gray 
scale  for  each  panel.  Each  panel  shows  a  reconstruction  area  of  5  mm  x  5  mm.  (a)  Cylinder,  radius 
2,5  mm,  with  an  internal  cylinder  of  radius  0.2  mm.  (b)  Cylinder,  radius  2.5  mm,  with  random 
internal  structure,  (c)  Tissue  structure,  with  variable  sound  speed  and  density,  from  chest  wall 


cross  section  5L  in  Ref.  56. 
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FIG.  7.  Three-dimensional  reconstructions  of  a  uniform  slab  with  contrast  7  =  0.01.  Each 


reconstruction  shows  an  isosurface  rendering  of  the  surface  7 m  =  0.0025.  Left:  single-frequency 
reconstruction.  Right:  time-domain  reconstruction. 
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Abstract 


A  finite-difference  time-domain  model  for  ultrasonic  pulse  propagation 
through  soft  tissue  has  been  extended  to  incorporate  absorption  effects  as 
well  as  longitudinal-wave  propagation  in  cartilage  and  bone.  This  extended 
model  has  been  used  to  simulate  ultrasonic  propagation  through  anatomically 
detailed  representations  of  chest  wall  structure.  The  inhomogeneous  chest  wall 
tissue  is  represented  by  two-dimensional  maps  determined  by  staining  chest 
wall  cross  sections  to  distinguish  between  tissue  types,  digitally  scanning  the 
stained  cross  sections,  and  mapping  each  pixel  of  the  scanned  images  to  fat, 
muscle,  connective  tissue,  cartilage,  or  bone.  Each  pixel  of  the  tissue  map 
is  then  assigned  a  sound  speed,  density,  and  absorption  value  determined 
from  published  measurements  and  assumed  to  be  representative  of  the  local 
tissue  type.  Computational  results  for  energy  level  fluctuations  and  arrival 
time  fluctuations  show  qualitative  agreement  with  measurements  performed 
on  the  same  specimens,  but  show  significantly  less  waveform  distortion  than 
measurements.  Visualization  of  simulated  tissue-ultrasound  interactions  in 
the  chest  wall  shows  possible  mechanisms  for  image  aberration  in  echocardio¬ 
graphy,  including  effects  associated  with  reflection  and  diffraction  caused  by 
rib  structures.  Comparison  of  distortion  effects  for  varying  pulse  center  fre¬ 
quencies  shows  that,  for  soft  tissue  paths  through  the  chest  wall,  energy  level 
and  waveform  distortion  increase  markedly  with  rising  ultrasonic  frequency 

and  that  arrival-time  fluctuations  increase  to  a  lesser  degree. 

43.80.Qf,  43.20.Fn,  43.58.Ta,  43.80.Cs 
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INTRODUCTION 


Echocardiography  is  widely  employed  for  diagnosis  of  cardiac  disease  including  valvular 
defects,  pericardial  elfusion,  and  wall  motion  abnormalities.1-3  Commonly,  echocardiography 
is  performed  noninvasively  through  the  chest  (transthoracic)  using  an  external  probe  placed 
on  the  chest  wall.  The  chest  wall,  however,  can  considerably  degrade  image  quality  because 
acoustic  paths  between  the  skin  and  heart  can  contain  ribs  and  cartilage  as  well  as  inhomo¬ 
geneous  muscle  and  fatty  tissue.  The  result  is  that  as  many  as  10-30%  of  patients  cannot 
be  successfully  imaged  with  present  transthoracic  techniques.4  This  limitation  of  transtho¬ 
racic  echocardiography  has  led  to  the  development  of  transesophageal  echocardiography,  in 
which  the  heart  is  imaged  by  a  transducer  inserted  within  the  esophagus.1-4  Although  trans¬ 
esophageal  echocardiography  provides  superior  image  quality,  resulting  in  high  diagnostic 
sensitivity  and  specificity,  the  invasiveness  of  the  procedure  is  accompanied  by  increased 
risk.3-6  For  this  reason,  improvement  in  the  noninvasive  transthoracic  approach  is  desirable, 
for  example,  by  the  development  of  methods  to  compensate  for  image  degradation  caused 
by  the  chest  wall. 

An  understanding  of  ultrasonic  aberration  produced  by  the  chest  wall  is  important  to 
the  development  of  appropriate  compensation  methods  for  transthoracic  ultrasonic  imaging. 
Direct  measurements  of  ultrasonic  distortion  produced  by  chest  wall  specimens7,8  have  been 
helpful.  The  study  reported  in  Ref.  7  shows  that  propagation  through  the  chest  wall  causes 
substantial  beam  distortion.  However,  that  study  did  not  distinguish  the  effect  of  soft  tissue 
from  effects  caused  by  rib  structures.  In  Ref.  8,  a  detailed  study  of  distortion  caused  by 
soft  tissue  paths  indicates  that  soft  tissue  distortion  in  the  chest  wall  is  substantially  less 
than  corresponding  distortion  in  the  human  abdominal  wall.  However,  distortion  caused 
by  ribs  was  only  treated  qualitatively  in  the  latter  study  because  the  physical  mechanisms 
of  rib-induced  distortion  could  not  be  adequately  described  by  the  method  reported  there. 
Although  a  model  of  ultrasound  propagation  in  the  chest  wall  has  previously  been  described,9 
that  model  is  based  on  coarse  depictions  of  chest  wall  morphology  including  homogeneous 
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tissue  layers  and  evenly-spaced,  uniformly-shaped  ribs.  These  previous  experiments  and 
simulations,  therefore,  have  left  gaps  in  knowledge  about  the  physical  causes  of  ultrasonic 
wavefront  distortion  caused  by  the  chest  wall. 

Recent  work  on  simulation  of  ultrasonic  pulse  propagation10-12  has  provided  insight 
about  the  wavefront  distortion  caused  by  the  human  abdominal  wall.  Although  these  studies 
have  provided  specific  information  about  the  relationships  between  soft  tissue  morphology 
and  ultrasonic  wavefront  distortion,  the  work  is  not  fully  applicable  to  distortion  caused  by 
the  human  chest  wall.  The  morphology  of  chest  wall  soft  tissue  has  been  hypothesized  to 
be  different  from  the  abdominal  wall  in  ways  that  affect  ultrasonic  wavefront  distortion.8 
Furthermore,  imaging  through  the  chest  wall  is  complicated  by  ribs  that  limit  the  usable 
acoustic  window  size  and  cause  scattering  and  reflection. 

The  study  reported  here  applies  quantitative  simulation  methods,  similar  to  those  pre¬ 
sented  in  Refs.  10  and  12,  to  anatomically  detailed  chest  wall  models  that  include  the  ribs.  To 
accurately  model  the  strong  losses  associated  with  propagation  through  bone  and  cartilage, 
the  finite-difference  method  described  in  Ref.  10  has  been  extended  to  include  absorption. 
Quantitative  descriptions  of  the  distortion  caused  by  soft  tissues  are  obtained  by  statistical 
analysis  of  simulated  distortion.  Visualizations  of  wavefronts  propagating  through  maps  of 
chest  cross  sections  provide  evidence  about  physical  relationships  between  wavefront  distor¬ 
tion  and  the  morphology  of  ribs  and  soft  tissue  structures  in  the  chest  wall.  Further  insight 
about  wavefront  distortion  mechanisms  is  provided  by  comparison  of  distortion  results  for 
incident  pulses  of  different  center  frequencies. 
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I.  THEORY 


Ultrasonic  pulse  propagation  through  the  human  chest  wall  is  modeled  here  using  the 
equations  of  motion  for  a  fluid  of  variable  sound  speed,  density,  and  absorption.  The  tissue 
is  assumed  motionless  except  for  small  acoustic  perturbations.  Absorption  is  included  using 
an  adaptation  of  the  Maxwell  solid  model,  in  which  all  absorption  effects  are  represented 
by  a  single  relaxation  time.  This  assumption  results  in  frequency-independent  absorption 
characteristics.13  For  such  a  fluid,  the  linearized  equations  of  mass  conservation,  momentum 
conservation,  and  state  can  be  combined  to  obtain  the  first-order,  two-dimensional,  coupled 
propagation  equations 

~  +  p(xi  v )  c(®»  y )2  v '  v(* >  y>  0  =  -<*(x,  y )  p{x,  y,  t ),  (l) 

p(x,  y)  -  +  Vp(g,  y,  t )  =  0.  (2) 


Here,  p(x,  y,  t)  is  the  acoustic  perturbation  in  fluid  pressure,  v(x,  y,  t)  is  the  vector  acoustic 
velocity,  p(x ,  y)  is  the  ambient  density,  c(x,  y)  is  the  ambient  sound  speed,  and  a(x,  y)  is  an 
absorption  coefficient  that  is  equivalent  to  the  inverse  of  a  spatially-dependent  relaxation 
time  r(x,y). 

The  absorption  coefficient  a,  defined  as  a  real  quantity,  is  related  to  the  energy  lost 
per  unit  length  as  follows.  The  propagation  equations  (1)  and  (2)  lead,  for  plane-wave 
propagation  of  the  form  p  =  e*(kx~ut\  to  the  dispersion  relation 

,  lo  L  ia  ,  ^ 

k  —  ~  t/l  4 - ,  (3) 

c  V  oj 

where  k  is  the  complex  wavenumber,  oj  is  the  (real)  radial  frequency  27t/,  and  c  is  the  (real) 
sound  speed.  The  imaginary  part  of  the  wavenumber  k  is  the  absorption  in  nepers  per  unit 
length.  Thus,  the  absorption  parameter  a  can  be  obtained  by  numerical  solution  of  the 
equation 


T _ r  /-i  loss  (dB/length)  T_ 

ImW=  20  log,,,  (e)  =I“ 


oj  L  ict 

-1/1 +  - 
c  V  oj  . 


(4) 
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Solution  of  Eq.  (4)  results  in  wavenumbers  having  a  real  part  that  differs  from  u ;/ c.  Since 
this  discrepancy  is  less  than  1%  over  the  range  of  tissue  properties  employed  in  the  present 
study,  use  of  absorption  coefficients  computed  from  Eq.  (4)  does  not  significantly  affect 
propagation  characteristics  except  by  adding  the  specified  absorption. 

Equations  (1)  and  (2)  were  solved  numerically  using  the  finite-difference  time-domain 
(FDTD)  method  described  in  Refs.  10  and  14.  This  method  is  a  two-step  MacCormack 
algorithm  that  is  fourth-order  accurate  in  space  and  second-order  accurate  in  time.  Further 
details  on  this  class  of  finite  difference  algorithms  can  be  found  in  Refs.  15-17. 

The  initial  condition  was  chosen  to  model  the  experimental  configuration  in  Ref.  8,  in 
which  a  spatially  broad,  nearly  planar  wavefront  was  emitted  from  a  wideband,  pulsed, 
unfocused  source  far  from  the  tissue  layer.  The  initial  wavefront  was  represented  in  the 
present  simulation  as  a  plane  wave  pulse  propagating  in  the  +y  direction: 


j>( X,  y,  0)  =  -  sin[&o(y  -  yb)] 
u(x,y,  0)  =  0,  and 
v(x,y,0)  =  — — — , 

r  c 


(5) 


where  the  wavenumber  k0  is  equal  to  27 r/0/c  for  a  center  frequency  of  /0  and  a  is  the  Gaussian 
parameter  of  the  pulse  temporal  envelope.  The  spatial  Gaussian  parameter  a  was  chosen 
to  simulate  the  bandwidth  of  the  pulse  used  in  the  experiments,  as  discussed  below  in  the 
Method. 

The  computational  configuration  is  analogous  to  that  described  in  Ref.  10.  The  domain 
of  computation  is  two-dimensional,  with  the  y  direction  taken  to  be  parallel  to  the  direction 
of  propagation  and  the  x  direction  parallel  to  the  initial  wavefront.  As  in  Ref.  10,  periodic 
boundary  conditions  were  applied  on  the  domain  edges  that  were  parallel  to  the  direction  of 
propagation,  while  radiation  boundary  conditions  were  applied  on  the  edges  perpendicular 
to  the  direction  of  propagation. 
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II.  METHOD 


This  study  employed  six  chest  wall  specimens  from  four  different  donors  between  79  and 
85  years  of  age.  One  specimen  (4L)  was  from  a  white  female,  while  the  others  were  from 
white  males.  All  specimens  were  obtained  from  autopsy,  stored  unfixed  at  —  20  °C,  and 
thawed  when  needed  for  study.  Wavefront  distortion  measurements  made  on  some  of  these 
specimens  as  well  as  other  specimens  are  described  in  Ref.  8.  The  nomenclature  employed 
here  for  the  cross  sections  corresponds  to  that  of  Ref.  8  for  the  whole  specimens  from  which 
the  cross  sections  were  taken;  each  cross  section  is  identified  by  a  donor  number  together 
with  “L”  or  “R”  to  indicate  whether  the  specimen  was  taken  from  the  left  or  right  side  of  the 
breastplate.  Additional  numbers  were  used  in  Ref.  8  to  indicate  the  intercostal  space  used 
in  each  measurement;  here,  letters  are  used  to  indicate  independent  acoustic  paths.  Four  of 
the  cross  sections  employed  here  (4L,  5L,  7L,  and  7R)  were  taken  from  specimens  employed 
in  Ref.  8.  Cross  sections  8L  and  8R  were  taken  from  specimens  through  which  ultrasonic 
transmission  was  measured  in  the  study  described  in  Ref.  8,  but  distortion  statistics  for 
these  specimens  were  not  reported  in  Ref.  8  because  of  limited  acoustic  windows. 

After  the  wavefront  distortion  measurements,  the  specimens  were  cut  into  cross  sections 
using  the  technique  described  in  Ref.  10.  The  cross  sections  were  fixed  and  stained  with  a 
modified  Gomori’s  trichrome  stain  according  to  the  procedure  detailed  in  Ref.  18,  so  that 
tissue  types  could  be  distinguished.  This  stain  colored  muscle  tissue  red  and  connective  tissue 
blue  while  leaving  the  fat  its  natural  color.  Calcified  tissue,  including  bone  and  cartilage 
in  the  current  specimens,  was  not  differentially  stained  by  this  technique,  but  the  natural 
contrast  between  bone,  cartilage,  and  marrow  was  sufficient  to  allow  tissue  mapping.  Full- 
color  300  d.p.i.  images  of  the  cross  sections  were  created  by  placing  each  stained  tissue  cross 
section  directly  onto  the  surface  of  a  flatbed  digital  scanner.  Image  editing  packages19  were 
used  to  manually  segment  the  cross  sectional  images,  i.e.,  to  map  the  images  into  regions 
that  corresponded  to  one  of  six  media.  The  media  were  water  (representing  water  external 
to  specimens  or  blood  inside  blood  vessels),  fat  (including  subcutaneous  fat,  fat  interlaced 
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within  muscle  layers,  and  marrow),  muscle,  connective  tissue  (including  skin,  septa,  and 
fasciae),  cartilage,  and  bone  (including  cortical  bone  and  trabeculae  within  cancellous  bone). 

The  six  segmented  tissue  maps  are  shown  in  Fig.  1.  All  of  the  cross  sections  contain  a 
layer  of  septated  subcutaneous  fat  below  the  skin.  Most  of  the  cross  sections  also  include  a 
layer  composed  primarily  of  the  major  pectoral  muscles  and  their  connective  fasciae  above 
the  ribs.  Between  the  ribs  are  regions  of  muscle  (internal  intercostal  and  external  intercostal 
groups)  interlaced  with  fat.  In  some  cases,  additional  thin  layers  of  fat  between  muscle  layers 
are  apparent.  Cross  sections  4L  and  7R  are  cut  along  the  intercostal  spaces  parallel  to  the 
ribs,  so  that  in  each  a  wide  cross  section  of  soft  tissue  appears.  Cross  sections  5L,  7L,  and 
8L  are  cut  perpendicular  to  the  ribs,  so  that  each  contains  soft-tissue  acoustic  paths  with 
width  equal  to  the  width  of  the  corresponding  intercostal  spaces.  Cross  section  8R  is  cut 
perpendicular  to  the  sternum  at  a  location  of  large  curvature  in  the  ribs,  so  that  the  ribs 
are  diagonally  sectioned.  Several  blood  vessels  appear  in  cross  sections  4L,  7L,  7R,  and  7R; 
the  largest  of  these  is  the  internal  mammary  artery. 

The  basic  structure  of  the  cross  sections  is  consistent  with  standard  descriptions  of  chest 
wall  anatomy.20,21  Ribs  appear  in  each  cross  section;  each  rib  is  composed  of  a  “costal 
cartilage”  near  the  sternum  (shown  in  most  of  the  cross  sections  considered  here)  attached 
to  a  “true  rib”  (composed  primarily  of  cancellous  bone)  at  the  edge  farther  from  the  sternum. 
In  the  cross  sections  considered  here,  the  costal  cartilages  are  primarily  composed  of  calcified 
cartilage,  surrounded  by  a  thin  layer  of  cortical  bone  (solid,  dense  bone  with  microscopic 
porous  structure),  which  in  turn  is  surrounded  by  the  periosteum,  a  thin  membrane  of 
connective  tissue.  Cross  sections  7L  and  7R  also  appear  to  contain  a  small  amount  of 
cortical  bone  in  the  central  portion  of  the  ribs.  This  phenomenon  may  be  associated  with 
advanced  calcification  known  to  occur  in  aging  humans.22  Cancellous  bone,  composed  of 
thin  trabeculae  that  forms  macroscopic  cells  filled  with  marrow,  is  seen  in  all  the  ribs  of 
cross  section  5L,  which  was  taken  at  a  distance  farther  from  the  sternum  so  that  the  true 
ribs,  rather  than  the  costal  cartilages,  were  included  in  this  cross  section.  Some  cancellous 
bone  is  also  apparent  within  portions  of  the  ribs  of  cross  sections  4L  and  8R.  In  each  case, 
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the  cancellous  bone  is  surrounded  by  a  thin  layer  of  cortical  bone  and  by  the  periosteum.  A 
portion  of  the  sternum,  composed  of  cancellous  bone  surrounded  by  cortical  bone,  is  visible 
at  the  left  side  of  cross  section  4L. 

The  density  and  sound  speed  arrays  needed  for  the  finite-difference  computation  were 
created  by  mapping  regions  of  the  segmented  tissue  images  to  reference  density  and  sound 
speed  values  for  the  five  tissue  types  and  water.  The  water  sound  speed  and  density  em¬ 
ployed  are  those  of  pure  water  at  body  temperature  (37.0 °C). 23,24  Sound  speeds  for  muscle 
and  fat  were  obtained  by  averaging  values  for  human  tissues  given  in  Refs.  25  and  26.  A 
representative  sound  speed  for  connective  tissue  was  determined  using  an  empirical  formula 
relating  collagen  content  to  ultrasonic  sound  speed27  together  with  a  measured  value  for  the 
collagen  content  of  human  skin.28  The  sound  speed  employed  for  bone  was  obtained  from  an 
average  of  values  reported  in  Ref.  29  for  longitudinal-wave  propagation  in  human  cortical 
bone.  The  sound  speed  used  here  for  cartilage  is  that  given  in  Ref.  30  as  quoted  in  Ref.  25. 
Density  values  for  soft  tissues  were  determined  from  Ref.  31  by  averaging  values  reported 
for  adipose  tissue,  skeletal  muscle,  and  skin,  respectively.  Density  values  employed  for  bone 
and  cartilage  are  average  values  from  Ref.  29. 

Attenuation  values  were  determined  from  measurements  summarized  in  Ref.  25  for  hu¬ 
man  fat  at  37°C,  human  bicep  muscle  at  37°  C,  human  skin  at  40°  C,  human  and  bovine 
cartilage  at  23°  C,  and  human  skull  (temperature  not  reported).  Attenuation  values  reported 
at  other  ultrasonic  frequencies  were  interpolated  (or,  for  the  skull  data,  extrapolated)  to  ob¬ 
tain  values  for  2.3  MHz  (corresponding  to  the  pulse  center  frequency  employed  here  and  in 
Ref.  8)  assuming  a  linear  dependence  of  attenuation  on  frequency.  The  attenuation  for  water 
was  estimated  by  extrapolating  frequency-  and  temperature-dependent  attenuation  values 
summarized  in  Ref.  32  to  2.3  MHz  and  37.0°  C.  The  values  of  tissue  parameters  employed 
in  the  present  study  are  given  in  Table  I. 

The  finite-difference  program  was  employed  to  compute  propagation  of  a  plane  wave  pulse 
through  each  scanned  cross  section  from  the  skin  to  the  peritoneal  membrane,  mimicking 
the  propagation  path  employed  in  the  distortion  measurements  of  Ref.  8.  The  spatial  step 
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size  of  the  finite-difference  grid  was  chosen  to  be  0.0442  mm,  or  1/15  wavelength  in  water 
at  the  center  frequency  of  2.3  MHz.  The  temporal  step  size  was  chosen  to  be  0.00725  /xs, 
for  an  optimal  Courant-Friedrichs-Levy  number  cAt/Ax  of  0.25. 17  The  Gaussian  parameter 
a  of  the  source  pulse  was  chosen  to  be  0.4766  mm  in  accordance  with  the  experimentally 
measured  pulse  bandwidth  (for  pulses  transmitted  through  a  water  path)  of  1.2  MHz.  Visual 
comparison  confirmed  that  the  simulated  pulse  closely  matched  the  measured  pulses  in  shape 
and  length. 

Each  simulation  was  performed  on  a  workstation  with  128  MB  of  random-access  memory. 
Finite-difference  grids  on  the  order  of  1500x1000  points  were  employed.  At  each  time  step, 
the  wave  field  was  updated  on  a  grid  subset  chosen  to  include  the  entire  support  of  the 
acoustic  wave  but  to  exclude  quiescent  regions.  The  entire  pressure  field  was  saved  as  a 
raster  image  at  intervals  of  0.725  /is  for  later  visualization.  The  computation  time  for  each 
simulation  was  on  the  order  of  five  hours.33 

Signals  were  recorded  for  8.62  fj, s  at  a  sampling  frequency  of  138  MHz  by  simulated 
apertures  with  dimensions  close  to  those  in  the  experimental  study  of  Ref.  8.  Positions 
of  all  simulated  apertures  employed  are  sketched  in  Fig.  1.  The  simulation  of  receiving 
elements  was  performed  by  integrating  the  locally-computed  pressure  over  the  element  pitch 
of  0.21  mm.  For  cross  sections  cut  parallel  to  the  ribs,  the  simulated  apertures  contained 
68  elements  for  an  aperture  width  of  14.28  mm.  For  cross  sections  cut  perpendicular  to 
the  ribs,  55  simulated  elements  were  used  to  form  11.55  mm  apertures.  Element  directivity 
effects  were  implicitly  included  by  the  integration  of  omnidirectional  sensitivity  functions 
over  the  width  of  each  element. 

A  one-dimensional  version  of  the  reference  waveform  method10,34  was  used  to  calculate 
the  arrival  time  of  the  pulse  at  each  receiving  position  in  the  simulation  data.  The  arrival 
time  fluctuations  across  the  receiving  aperture  caused  by  each  cross  section  were  calculated 
by  subtracting  a  linear  fit  from  these  calculated  arrival  times.  Energy  level  fluctuations 
in  the  data  were  calculated  by  summing  the  squared  amplitudes  of  each  waveform  over  a 
2.4  n s  window  that  isolated  the  main  pulse,  converting  to  decibel  units,  and  subtracting 
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the  best  linear  fit  from  the  resulting  values.  As  for  polynomial  fits  previously  employed  in 
wavefront  distortion  measurements,8  the  purpose  of  the  linear  fit  removal  in  each  case  was  to 
compensate  for  gross  changes  in  tissue  thickness  across  the  array.  Variations  in  pulse  shape 
across  the  aperture  were  evaluated  using  the  waveform  similarity  factor;34  this  quantity, 
which  can  be  considered  a  generalized  cross-correlation  coefficient,  has  a  maximum  of  unity 
when  all  received  waveforms  are  identically  shaped. 

To  test  the  frequency  dependence  of  chest  wall  wavefront  distortion,  propagation  through 
eight  portions  of  specimens,  each  containing  only  soft  tissue,  was  also  computed  for  wave- 
fronts  having  center  frequencies  of  1.6  and  3.0  MHz.  In  each  case,  the  initial  wavefront  was 
chosen  to  have  the  same  temporal  envelope  as  above.  The  absorption  coefficient  at  these 
frequencies  for  each  tissue  type  was  extrapolated  from  the  value  employed  at  2.3  MHz  using 
the  assumption  that  absorption  depended  linearly  on  the  center  frequency.  The  spatial  and 
temporal  sampling  rates  were  also  varied  in  inverse  proportion  to  the  pulse  center  frequency. 
All  runs  were  otherwise  identical  in  configuration  and  processing  to  those  described  above. 
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III.  RESULTS 


Wavefront  distortion  results  for  13  soft  tissue  paths  (*.e.,  paths  in  which  wavefront  dis¬ 
tortion  was  not  significantly  influenced  by  the  ribs)  are  shown  in  Table  II.  These  results 
indicate  that  soft  tissue  paths  cause  a  wide  range  of  wavefront  distortion  effects  depending 
on  the  specific  morphology  of  each  path.  For  instance,  path  7R-c  causes  arrival  time  and 
energy  level  fluctuations  that  are  more  than  twice  the  magnitude  of  those  caused  by  the 
adjacent  path  7R-d.  This  difference  is  thought  to  arise  from  morphological  features,  includ¬ 
ing  muscle  tissue  with  interlaced  fat  and  a  large  amount  of  connective  tissue,  of  the  tissue 
within  path  7R-c.  Also  notable  is  that  the  specimen  thickness  does  not  closely  correspond 
to  variations  in  distortion.  The  largest  rms  arrival  time  fluctuation  and  lowest  waveform 
similarity  factor,  for  example,  are  caused  by  path  4L-c,  which  has  an  average  thickness  less 
than  the  mean  for  all  the  tissue  paths. 

Wavefront  distortion  statistics  for  the  13  soft  tissue  paths  are  graphically  summarized  in 
Fig.  2  together  with  corresponding  statistics  for  all  of  the  soft  tissue  measurements  reported 
in  Ref.  8.  This  comparison  indicates  that  wavefront  distortion  caused  by  soft  tissues  in  the 
chest  wall  simulations  is  comparable  to  measured  distortion.  Arrival  time  fluctuations  and 
energy  level  fluctuations  for  simulated  distortion  are  slightly  less  than  measured  values,  but 
mean  values  of  both  fluctuations  for  the  simulations  fall  well  within  one  standard  deviation  of 
the  corresponding  mean  fluctuation  for  the  measurements.  The  waveform  similarity  factor, 
however,  is  substantially  higher  for  simulations  than  measurements,  indicating  that  sim¬ 
ulated  waveforms  were  distorted  considerably  less  than  measured  waveforms.  Correlation 
lengths  for  the  simulated  distortions  are  somewhat  less  than  measured  values.  However, 
consistent  with  measurements,  the  mean  correlation  length  of  the  simulated  arrival  time 
fluctuations  is  greater  than  that  for  the  simulated  energy  level  fluctuations. 

As  in  Ref.  8,  rib  structures  were  found  to  cause  much  more  distortion  than  soft  tissue 
alone.  The  varied  nature  of  distortion  caused  by  rib  effects  is  illustrated  in  Fig.  3,  which 
shows  three  representative  sets  of  measured  signals  for  specimen  8L.  These  measurements 
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were  obtained  by  the  methods  reported  in  Ref.  8.  The  first  panel  shows  96  adjacent  measured 
signals,  along  the  array  direction  (approximately  parallel  to  the  ribs)  for  propagation  through 
a  tissue  path  within  an  intercostal  space.  The  signals  are  not  severely  distorted;  secondary 
arrivals  are  discernible,  but  are  of  lower  amplitude  than  the  main  arrival.  The  second  panel 
shows  96  measured  signals  for  an  elevation  over  a  rib.  Here,  all  signals  are  severely  distorted. 
Multiple  arrivals,  as  well  as  high- amplitude  spatially-random  fluctuations,  are  seen.  The 
third  panel  shows  50  measured  signals  along  the  elevation  direction  (perpendicular  to  the 
ribs),  centered  over  the  soft  tissue  between  the  ribs.  Here,  the  main  wavefront  is  curved 
rather  than  straight,  an  additional  arrival  behind  the  main  wavefront  is  seen,  and  portions 
of  the  signals  from  over  the  ribs  (at  both  edges  of  the  panel)  are  advanced  relative  to  the 
signals  from  the  central  soft  tissue  region. 

The  present  simulations  allow  more  detailed  qualitative  and  quantitative  investigation 
of  rib  effects  than  were  possible  from  the  previous  measurements.  Propagation  through  two 
rib-influenced  paths  is  illustrated  in  Figs.  4  and  5  by  images  shown  in  formats  similar  to 
visualizations  of  propagation  through  soft  tissue  in  Ref.  10.  Simulated  ultrasonic  pulses  are 
superimposed  on  portions  of  the  tissue  maps  from  Fig.  1. 

Figure  4  shows  propagation  through  a  thin  rib,  composed  chiefly  of  cancellous  bone,  in 
cross  section  5L  (corresponding  approximately  to  path  5L-b).  A  strong  reflection  occurs  at 
the  first  interface  between  bone  and  soft  tissue,  removing  a  substantial  amount  of  energy 
from  the  main  wavefront.  The  small,  high-contrast  trabeculae  within  the  rib  cause  consid¬ 
erable  scattering,  as  can  be  observed  in  panel  (b)  of  Fig.  4.  The  scattering  causes  random 
fluctuations  behind  the  main  wavefront;  these  fluctuations  somewhat  resemble  those  seen  in 
the  measured  data  of  Fig.  3(b).  After  passing  through  the  rib,  as  seen  in  panels  (c)  and  (d) 
of  Fig.  4,  the  central  portion  of  the  wavefront  shows  substantial  attenuation  and  distortion. 
However,  the  average  arrival  time  of  the  wavefront  is  not  greatly  changed  by  propagation 
through  the  rib,  but  is  advanced  by  only  about  one-half  period.  This  phenomenon  appar¬ 
ently  occurs  because  the  influence  of  the  “slow”  marrow  (modeled  here  as  fat)  counteracts 
the  influence  of  the  “fast”  trabeculae.  Noteworthy  is  that  the  predominant  ultrasonic  wave- 
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length  has  increased  after  propagation  through  the  rib,  so  that  the  effective  center  frequency 
of  the  wavefront  has  been  lowered.  Since  the  absorption  model  used  in  the  present  study 
includes  only  frequency-independent  absorption,  the  loss  of  short-wavelength  components  in 
this  simulation  results  only  from  frequency-dependent  scattering  caused  by  the  trabeculae. 

Propagation  within  path  8L-b,  which  includes  two  larger  ribs  and  the  corresponding 
intercostal  space,  is  illustrated  in  Figure  5.  At  the  position  of  the  cross  section,  these 
ribs  are  composed  primarily  of  cartilage  and  surrounded  by  a  thin  layer  of  cortical  bone. 
Since  the  cartilage  and  bone  of  these  ribs  is  modeled  as  homogeneous,  small-scale  scattering 
within  these  tissues  did  not  occur  in  this  simulation.  Instead,  the  wavefront  is  reflected  from 
interfaces  between  cartilage,  bone,  and  soft  tissue. 

The  visualization  shown  in  Fig.  5  provides  physical  reasons  for  all  the  rib-related  distor¬ 
tion  phenomena  seen  in  the  measured  data  of  Fig.  3(c).  The  wavefronts  propagating  through 
the  ribs  show  greater  attenuation  than  that  in  Fig.  4,  both  because  of  the  high  absorption 
of  the  ribs  and  because  of  the  reflections  noted  above.  These  wavefronts  are  also  advanced 
relative  to  the  wavefront  propagating  through  the  intercostal  space,  because  of  the  higher 
sound  speed  of  both  bone  and  cartilage.  The  wavefront  propagating  through  the  intercostal 
space  is  distorted  somewhat  by  the  inhomogeneous  soft  tissue  path,  as  can  be  observed  in 
panels  (b)  and  (c).  However,  much  greater  distortion  results  from  interaction  between  the 
wavefront  and  the  ribs.  A  rightward-propagating  reflection,  seen  in  panels  (b)  and  (c),  com¬ 
bines  with  the  main  wavefront  in  panel  (d)  to  result  in  severe  distortion  at  the  right  side  of 
the  central  wavefront.  Furthermore,  diffraction  from  the  edges  of  the  ribs  results  in  large 
curvature  of  the  soft-tissue  wavefront. 

Distortion  and  attenuation  statistics  for  a  variety  of  rib-influenced  paths  are  shown  in 
Table  III.  Footnotes  in  Table  III  indicate  physical  causes  of  distortion  present  within  each 
path.  A  variety  of  distortion  and  attenuation  mechanisms  are  illustrated.  Propagation 
through  small  intercostal  spaces  (paths  4L-a,  8L-b,  8L-f,  and  7R-a)  causes  diffraction  effects 
that  introduce  substantial  curvature  into  the  wavefront,  as  seen  in  Fig.  5.  This  large-scale 
wavefront  curvature  is  associated  with  large  arrival  time  fluctuation  values  although  the 
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wavefronts  generally  appear  to  be  locally  smooth.  Interference  between  directly-transmitted 
and  rib-reflected  wavefronts  (paths  4L-a,  8L-b,  8L-d,  8L-f,  and  7R-a)  introduces  arrival 
time,  energy  level,  and  waveform  distortion  substantially  greater  than  that  for  soft  tissue 
paths  without  ribs.  Propagation  through  cancellous  bone  (paths  4L-a,  4L-b,  5L-b,  and 
8R-c)  results  in  considerable  attenuation  and  large  waveform  distortion,  while  propagation 
through  cortical  bone  and  cartilage  (paths  4L-a,  4L-b,  8L-a,  8L-c,  8L-e,  8L-g,  7L-c,  7R-a, 
7R-b,  and  8R-c)  results  in  even  larger  attenuation  but  smaller  distortion.  Where  bone  is 
embedded  within  cartilage  (paths  7-Lc  and  7R-b),  additional  scattering  also  occurs.  For  the 
path  including  a  large  bone  inclusion  (paths  7R-b),  this  scattering  results  in  extremely  high 
energy  level  and  waveform  distortion. 

Frequency-dependent  wavefront  distortion  statistics  are  summarized  in  Fig.  6.  Tissue 
paths  used  for  these  computations,  none  of  which  include  rib  structures,  are  those  labeled 
4L-d,  4L-f,  5L-a,  5L-c,  8R-a,  8R-b,  7L-a,  and  7L-b  in  Fig.  1.  The  results  shown  in  Fig.  6 
indicate  that  arrival  time  fluctuations,  energy  level  fluctuations,  and  waveform  distortion 
all  become  more  severe  with  increasing  pulse  frequency.  The  most  dramatic  change  is  in 
the  energy  level  distortion;  on  average,  the  rms  energy  level  fluctuations  for  the  3.0  MHz 
signals  are  2.3  times  those  for  the  1.6  MHz  signals.  Correlation  lengths  of  both  arrival 
time  and  energy  level  fluctuations  decrease  with  frequency,  so  that  the  predominant  length 
scales  of  ultrasonic  wavefront  distortion  are  seen  to  decrease  with  the  ultrasonic  wavelength. 
As  with  the  rms  distortion  statistics,  the  most  dramatic  frequency-dependent  change  is  in 
the  energy  level  fluctuations.  Still,  even  the  high-frequency  pulses  here  show  substantially 
smaller  distortion  than  that  previously  observed  in  experiments  and  simulations  for  the 
human  abdominal  wall. 10-12,35 
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IV.  DISCUSSION 


As  with  earlier  simulations  of  propagation  through  tissue,10,12  the  current  study  shows 
qualitative  agreement  with  measured  wavefront  distortion  results  for  similar  specimens.8 
However,  the  accuracy  of  the  present  model  is  limited  by  simplifications  of  true  tissue  struc¬ 
ture.  In  particular,  the  computational  model  here  does  not  account  for  property  variations 
within  tissue  types,  tissue  microstructure,  or  three-dimensional  tissue  structure.  These  lim¬ 
itations  are  discussed,  with  respect  to  soft  tissues,  in  Ref.  10. 

The  modeling  of  ribs  adds  additional  complication.  In  the  current  study,  individual 
trabeculae  were  assumed  to  be  composed  of  tissue  having  properties  identical  to  cortical 
bone,  an  assumption  known  as  Wolff’s  hypothesis.36  The  validity  of  this  hypothesis  has  been 
questioned;37,38  however,  measured  elastic  properties  of  individual  trabeculae  vary  widely37,38 
and  recent  work39  has  provided  support  for  Wolff’s  hypothesis.  Thus,  the  properties  em¬ 
ployed  here  for  trabecular  bone  can  be  regarded  as  a  reasonable  order-of-magnitude  estimate. 
Likewise,  the  modeling  of  marrow  as  fat  tissue  is  a  simplifying  assumption  that  may  have 
limited  validity,  although  available  data  suggest  that  the  density  and  sound  speed  of  marrow 
are  close  to  those  for  other  adipose  tissues.29  In  addition,  the  present  model  for  cartilage 
is  based  on  measurements  of  normal  cartilage,  while  the  cartilage  present  in  the  specimens 
employed  here  was  calcified  due  to  the  age  of  the  donors.  However,  density  measurements 
made  on  eight  representative  samples  of  calcified  cartilage  (two  from  specimen  7R,  four 
from  specimen  1R,8  and  two  from  an  unused  specimen)  resulted  in  an  average  density  of 
0.00111  kg/m3,  which  is  different  by  only  1%  from  the  density  assumed  here.  Since  sound 
speed  in  calcified  tissue  has  been  empirically  shown  to  be  directly  related  to  density,40,41  this 
small  change  in  density  suggests  that  the  acoustic  properties  of  the  calcified  cartilage  in  our 
specimens  is  close  to  that  for  normal  cartilage. 

The  computations  reported  here  model  the  chest  wall  as  a  fluid  of  variable  sound  speed, 
density,  and  compressibility.  This  model  implicitly  neglects  shear  wave  propagation.  The 
neglect  of  shear  waves  in  soft  tissues  is  believed  to  be  justified  because  absorption  of  shear 
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waves  in  soft  tissues  is  much  greater  than  absorption  of  longitudinal  waves.42,43  In  calcified 
tissues,  however,  significant  shear  waves  are  known  to  be  generated.44,45  In  the  current 
scattering  configuration,  some  shear  waves  are  likely  generated  wherever  the  rib  surface  is 
far  from  parallel  to  the  wavefront.  However,  since  shear  wave  absorption  has  been  found  to 
be  somewhat  larger  than  longitudinal  wave  absorption  for  ultrasonic  propagation  in  bone,44 
the  significance  of  shear-wave  propagation  within  bone  on  transmitted  ultrasonic  wavefronts 
is  questionable.  For  this  reason,  omission  of  nonlongitudinal  waves  in  the  present  study, 
as  in  another  computational  study  of  ultrasonic  scattering  from  bone,46  is  believed  to  be 
justified;  however,  further  study  would  be  required  to  confirm  this  assumption. 

The  absence  of  frequency-dependent  absorption  is  a  possible  source  of  error  in  the  present 
estimates  of  total  tissue  attenuation,  energy  level  fluctuations,  and  waveform  distortion. 
However,  since  absorption  in  tissue  increases  approximately  linearly  with  frequency,  lower 
absorption  for  frequency  components  below  the  pulse  center  frequency  would  nearly  cancel 
higher  absorption  for  frequency  components  above  the  center  frequency,  so  that  the  average 
absorption  incurred  by  a  wideband  pulse  should  still  be  computed  with  fair  accuracy.  For 
this  reason,  the  absence  of  frequency-dependent  absorption  in  the  calculations  reported  here 
is  not  considered  to  be  a  significant  source  of  error  in  the  computed  attenuation  or  energy 
level  fluctuation  curves.  Still,  inclusion  of  frequency-dependent  absorption  would  result  in 
additional  waveform  distortion  effects.  The  lack  of  this  effect  is  a  likely  reason  for  the 
lower  waveform  distortion  (higher  waveform  similarity  factors)  obtained  from  simulations 
as  compared  to  measurements.  However,  the  absence  of  frequency-dependent  absorption 
effects  allowed  frequency-dependent  scattering  effects  to  be  clearly  quantified  separately 
from  absorption  effects. 

Although  the  simulations  were  planned  to  match  the  measurements  of  Ref.  8  closely, 
a  number  of  differences  remain.  The  most  important  of  these,  as  discussed  in  Ref.  10,  is 
that  the  simulations  were  performed  using  a  two-dimensional  tissue  model  while  the  mea¬ 
surements  were  inherently  three-dimensional.  Other  differences  include  details  of  the  source 
waveform  and  wavefront  shape,  variations  in  the  specimen  orientations  and  the  regions  in- 
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terrogated,  and  variations  in  the  distance  between  the  specimen  and  the  real  or  simulated 
receiving  aperture.  All  of  these  differences  could  contribute  to  discrepancies  between  mea¬ 
surements  and  simulations. 

In  general,  most  of  the  simplifying  assumptions  in  the  present  tissue  model  are  likely  to 
result  in  underestimation  of  wavefront  distortion  produced  by  the  human  chest  wall.  In¬ 
corporation  of  tissue  microstructure,  spatially-dependent  acoustic  properties  for  each  tissue 
type,  shear  wave  propagation  in  bone  and  cartilage,  three-dimensional  propagation,  and 
frequency-dependent  absorption  could  all  result  in  greater  spatial  and  temporal  variations 
in  the  propagating  acoustic  fields,  so  that  these  features  could  produce  simulated  distortion 
with  characteristics  closer  to  measurements.  For  this  reason,  distortion  statistics  computed 
using  the  present  tissue  model  should  be  interpreted  as  lower  limits  for  the  statistics  of 
distortion  occurring  in  real  chest  wall  tissue. 

Additionally,  some  of  the  discrepancy  between  simulated  and  measured  distortion  may 
be  explained  by  the  nonuniform  characteristics  of  the  receiving  transducer  employed  in  the 
measurements.8  The  water-path  measurements  reported  in  Ref.  8  show  arrival  time  fluctu¬ 
ations  (mean  2.21  ns)  and  energy  level  fluctuations  (mean  0.36  dB);  although  small,  these 
fluctuations  are  comparable  to  the  difference  between  the  average  measured  and  simulated 
fluctuations.  Thus,  compensation  for  arrival  time  and  energy  level  fluctuations  due  to  trans¬ 
ducer  irregularities  could  reduce  measured  distortion  to  levels  closer  to  the  simulations. 
Also,  the  waveform  similarity  factor  for  water  path  measurements  was  0.991, 8  which  indi¬ 
cates  greater  waveform  distortion  than  the  average  value  of  0.995  computed  here  for  soft 
tissue  paths.  Thus,  compensation  of  the  measured  data  for  transducer  impulse-response  vari¬ 
ations  could  raise  the  measured  waveform  similarity  factor  to  a  value  in  closer  agreement 
with  simulations. 

Previous  experimental  measurements  of  wavefront  distortion  caused  by  the  human  chest 
wall8  have  suggested  that  distortion  caused  by  chest  wall  soft  tissues  is  less  severe  than  that 
caused  by  the  human  abdominal  wall.11,35  This  difference  has  been  observed  to  occur  even 
though  average  specimen  thicknesses  were  comparable  in  chest  wall8  and  abdominal  wall11,35 
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measurements.  The  present  study  provides  support  for  these  results;  arrival  time  and  energy 
level  distortion  by  the  chest  wall  was  found  here  to  be  smaller  than  that  produced  by  the 
abdominal  wall  in  previous  simulation  studies.10,12  For  the  simulations,  this  difference  may 
be  partially  explained  by  the  fact  that  the  chest  wall  specimens  employed  here  are  thinner 
on  average  (mean  thickness  17.5  mm)  than  the  abdominal  wall  cross  sections  employed  in 
Refs.  10  and  Ref.  12  (mean  thickness  26.7  mm).  Another  possible  partial  explanation  is 
that  the  pulse  center  frequency  employed  in  abdominal  wall  measurements  and  simulations 
was  3.75  MHz,  significantly  higher  than  the  center  frequency  of  2.3  MHz  for  the  chest  wall 
measurements  and  simulations.  Differences  in  pulse  frequency  and  specimen  thickness  may 
explain  the  discrepancy  in  energy  level  distortion  between  the  abdominal  wall  and  chest  wall, 
but  do  not  fully  explain  the  discrepancy  in  arrival  time  distortion  results.  For  instance, 
the  mean  arrival  time  and  energy  level  fluctuations  per  unit  length  are  1.02  ns/mm  and 
0.083  dB/mm  for  the  present  study  vs.  1.96  ns/mm  and  0.105  dB/mm  for  the  abdominal 
wall  cross  sections  of  Ref.  10  and  12.  Arrival  time  distortion  was  shown  here  to  increase  only 
subtly  with  increasing  pulse  frequency,  so  that  this  discrepancy  in  arrival  time  fluctuations  is 
not  fully  explained  by  pulse  frequency  differences.  However,  energy  level  fluctuations  increase 
markedly  with  frequency  for  chest  wall  tissue.  Thus,  for  equal  ultrasonic  pulse  frequencies, 
chest  wall  tissue  should  cause  energy  level  distortion  per  unit  length  comparable  to  that 
caused  by  abdominal  wall  tissue. 

It  was  suggested  in  Ref.  8  that  chest  wall  morphology  may  differ  from  abdominal  morphol¬ 
ogy  in  a  manner  that  results  in  smaller  ultrasonic  wavefront  distortion.  The  cross  sections 
employed  here  can  be  compared  with  those  employed  in  Ref.  10  and  Ref.  12  to  evaluate 
the  importance  of  morphological  differences  between  chest  wall  and  abdominal  wall  tissue. 
One  difference  between  the  two  groups  of  cross  sections  is  the  nature  of  the  subcutaneous 
fat  layers.  The  abdominal  wall  cross  sections  generally  contain  thicker  fat  layers,  containing 
many  more  lobular  structures  than  the  chest  wall  cross  sections.  Since  the  high  contrast 
between  septa  and  fat  causes  substantial  ultrasonic  scattering,10-12  this  morphological  dif¬ 
ference  is  likely  to  result  in  lower  overall  energy  level  and  waveform  distortion  for  chest 
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wall  tissue  (although,  as  discussed  above,  the  energy  level  distortion  per  unit  propagation 
length  should  be  comparable).  Also,  the  abdominal  wall  and  chest  wall  cross  sections  have 
markedly  different  structure  within  the  muscle  layers  that  occur  below  the  subcutaneous 
fat.  The  abdominal  wall  cross  sections  have  many  large-scale  features  due  to  aponeuroses 
(interfaces  between  muscle  groups,  composed  of  connective  tissue  and  fat)  and  large  fatty 
regions.  These  large-scale  features  cause  large  wavefront  fluctuations  that  are  associated 
with  large  rms  arrival  time  fluctuations.10,12  In  contrast,  muscle  layers  of  the  chest  wall  cross 
sections  considered  here  contain  primarily  smaller-scale  structure  associated  with  small  is¬ 
lands  of  interlaced  fatty  tissue.  This  morphological  difference  may  result  in  lower  large-scale 
arrival  time  fluctuations  but  significant  energy  level  fluctuations  associated  with  scattering, 
consistent  with  the  differences  between  distortion  caused  by  soft  tissues  in  the  abdominal 
wall  and  the  chest  wall. 

The  present  results  for  frequency  dependence  of  distortion  provide  further  insight  into 
the  importance  of  scattering  effects  relative  to  large-scale  structure  in  wavefront  distortion 
caused  by  soft  tissues.  If  wavefront  distortion  in  the  chest  wall  were  caused  only  by  large- 
scale  tissue  structures,  wavefront  distortion  to  be  roughly  independent  of  frequency,  since 
propagation  effects  are  independent  of  frequency  in  the  geometric  acoustics  limit.  How¬ 
ever,  distortion  caused  by  scattering  effects  should  increase  with  the  pulse  frequency  for 
inhomogeneities  of  size  comparable  to  the  wavelength.  Previous  simulation  and  experimen¬ 
tal  studies10-12  on  distortion  caused  by  the  human  abdominal  wall  have  suggested  that 
energy  level  fluctuations  and  waveform  distortion  are  generally  associated  with  scattering 
effects,  while  arrival  time  fluctuations  are  predominantly  caused  by  large-scale  path  length 
differences.  The  present  results,  while  consistent  with  those  conclusions,  indicate  that  scat¬ 
tering  plays  a  role  in  all  types  of  distortion  considered  here.  Since  energy  level  fluctuations 
and  waveform  similarity  factors  exhibit  more  dramatic  increases  in  distortion  with  increas¬ 
ing  pulse  frequency,  the  present  results  suggest  that  scattering  is  of  primary  importance  in 
causing  energy  level  and  waveform  distortion  and  of  secondary  importance  in  causing  arrival 
time  distortion. 
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These  results  can  be  employed  to  evaluate  the  potential  of  various  approaches  to  im¬ 
prove  echocardiographic  imaging.  Available  acoustic  windows  for  transthoracic  imaging  are 
severely  limited  by  the  presence  of  the  ribs,  so  that  image  quality  cannot  be  significantly 
improved  by  increase  of  aperture  size.  The  present  results  also  indicate  that  use  of  higher- 
frequency  probes  may  provide  less  benefit  than  expected  because  of  frequency-dependent 
scattering  in  the  chest  wall. 

For  these  reasons,  aberration  correction  methods  are  potentially  important  in  transtho¬ 
racic  echocardiography,  particularly  for  higher-frequency  imaging.  The  frequency-dependent 
distortion  results  reported  here  suggest  that  distortion  models  employing  single  phase  screens 
may  be  of  some  benefit  for  aberration  correction  in  echocardiography  through  soft  tissue 
paths.  The  relatively  weak  dependence  of  arrival  time  fluctuations  on  pulse  frequency  sug¬ 
gests  that  a  large  portion  of  arrival  time  variations  are  caused  by  tissue  structures  too  large 
to  cause  significant  frequency-dependent  scattering  effects.  Similar  conclusions  regarding 
the  importance  of  large-scale  structure  to  arrival  time  fluctuations  have  also  been  drawn 
from  results  presented  in  Refs.  10  and  12. 

Still,  the  present  results,  like  those  from  earlier  studies,10-12  suggest  that  single  phase 
screens  will  not  provide  complete  correction  for  distortion  caused  by  soft  tissues.  In  partic¬ 
ular,  methods  employing  single  phase  screens  will  not  completely  remove  distortion  caused 
by  scattering.  The  sharp  increase  of  amplitude  and  waveform  distortion  with  frequency,  as 
well  as  the  moderate  increase  of  arrival  time  distortion  with  frequency,  indicate  that  scat¬ 
tering  effects  become  much  more  important  to  ultrasonic  aberration  as  imaging  frequencies 
increase.  Furthermore,  phase  screen  models  do  not  inherently  account  for  distortion  caused 
by  rib  structures,  shown  here  to  produce  diffraction,  reflection,  and  scattering.  Thus,  any 
attempted  correction  using  only  phase  screen  models  is  likely  to  provide  little  improvement 
in  the  presence  of  strong  rib-induced  effects. 

Other  correction  models  that  incorporate  rib  structures  may  provide  greater  image  im¬ 
provements  for  the  distortions  most  important  to  echocardiography.  Processing  wavefronts 
with  techniques  such  as  angular  spectrum  filtering  can  remove  some  spurious  arrivals,47  al- 
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though  such  computations  may  be  difficult  to  incorporate  into  a  general  correction  algorithm. 
Other  possible  methods  include  those  incorporating  models  of  tissue  structure.  Models  in¬ 
corporating  ray  acoustics9  may  provide  improvement,  but  implicitly  neglect  diffraction  and 
scattering  effects,  so  that  aberration  correction  would  be  incomplete,  particularly  for  small 
intercostal  spaces.  A  more  complete  aberration  correction  method  could  employ  synthetic 
focusing  using  full-wave  numerical  computation  of  acoustic  fields  within  sufficiently  accu¬ 
rate  models  of  tissue  structure.  This  method  has  been  implemented,  within  the  context  of  a 
quantitative  frequency-domain  inverse  scattering  method,  in  Ref.  48.  However,  the  results 
presented  here  indicate  that  distortion  caused  by  soft  tissue  and  rib  structures  varies  widely 
based  on  morphological  variations  between  (and  within)  individuals.  Thus,  for  any  general 
correction  method  employing  models  of  tissue  structure,  separate  models  of  tissue  structure 
must  be  constructed  for  each  region  of  interest. 
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V.  CONCLUSIONS 


A  computational  study  of  ultrasonic  propagation  through  the  chest  wall,  including  tissue- 
dependent  absorption  as  well  as  detailed  anatomical  cross  sections,  has  been  presented.  For 
soft-tissue  paths,  computational  results  for  arrival  time  distortion,  energy  level  distortion, 
and  correlation  lengths  of  these  distortions  are  comparable  to  those  reported  in  previous 
chest  wall  measurements.  Both  simulations  and  measurements  indicate  that  arrival  time 
distortion  and  energy  level  distortion  caused  by  soft  tissues  in  the  human  chest  wall  is 
smaller  than  that  caused  by  the  human  abdominal  wall.  Differences  in  morphology  between 
the  abdominal  wall  and  the  chest  wall  provide  a  probable  explanation  for  this  difference. 

Distortion  caused  by  rib  structures  is  much  more  severe  than  that  caused  by  soft  tissues. 
Reflections  and  diffraction  from  rib  structures  complicate  wavefronts  that  travel  through 
soft  tissue  paths  adjacent  to  ribs  and  can  cause  arrival  time  and  energy  level  fluctuations 
much  greater  than  those  induced  by  soft  tissue  structures.  Wavefronts  propagating  directly 
through  rib  structures  are  attenuated  by  both  internal  absorption  and  reflection  at  interfaces 
between  bone,  cartilage,  and  soft  tissue.  Internal  scattering  within  rib  structures  causes 
distortion  phenomena  that  include  severe  waveform  and  energy  level  distortion,  additional 
attenuation,  and  lowering  of  the  effective  frequency  for  the  transmitted  pulse.  The  strong 
dependence  of  distortion  on  the  morphological  details  of  rib  structures  presents  a  major 
challenge  for  aberration  correction  in  echocardiography. 

Simulation  of  propagation  through  soft-tissue  paths  using  three  different  pulse  frequencies 
has  indicated  that  the  distortion  types  investigated  here  have  different  frequency  dependence. 
Arrival  time  fluctuations  increase  subtly  with  frequency,  while  energy  level  and  waveform 
distortion  increase  greatly.  Thus,  a  substantial  portion  of  arrival  time  fluctuations  in  the 
chest  wall  may  be  explained  by  large-scale  tissue  variations,  but  some  arrival  time  distortion 
and  most  energy  level  and  waveform  distortion  apparently  result  from  scattering.  Thus, 
correction  of  wavefront  distortion  caused  by  soft  tissues  should  become  both  more  important 
and  more  challenging  as  pulse  frequencies  employed  in  imaging  systems  are  increased. 
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TABLES 


TABLE  I.  Assumed  physical  properties  for  each  tissue  type  employed  in  the  simulations. 


Tissue 

Type 

Sound  Speed 
(mm/fj,s) 

Density 

(g/cc) 

Absorption 

(dB/mm) 

Water 

1.524 

0.993 

0.0007 

Fat 

1.478 

0.950 

0.12 

Muscle 

1.547 

1.050 

0.21 

Connective 

1.613 

1.120 

0.37 

Cartilage 

1.665 

1.098 

0.97 

Bone 

3.540 

1.990 

4.37 

TABLE  II.  Statistics  of  wavefront  distortion  caused  by  thirteen  soft-tissue  paths  within  chest 
wall  cross  sections.  The  statistics  shown  include  the  average  specimen  thickness  for  the  tissue  path 
considered,  rms  values  and  correlation  lengths  (CL)  of  the  arrival  time  fluctuations  (ATF)  and  the 
energy  level  fluctuations  (ELF),  the  waveform  similarity  factor  (WSF),  and  the  total  attenuation. 


Path 

Thickness 

(mm) 

ATF 

rms 

(ns) 

CL 

(mm) 

ELF 

rms 

(dB) 

CL 

(mm) 

WSF 

Attenuation 

(dB) 

4L-c 

15.4 

32.0 

0.60 

1.98 

1.68 

0.981 

5.62 

4L-d 

12.7 

10.0 

2.58 

0.46 

1.23 

0.999 

4.08 

4L-e 

16.0 

10.0 

1.37 

1.61 

1.74 

0.998 

5.26 

4L-f 

17.0 

17.3 

2.48 

0.92 

1.61 

0.999 

5.33 

5L-a 

11.0 

11.6 

0.95 

1.51 

1.13 

0.991 

4.29 

5L-c 

15.0 

14.8 

1.03 

1.15 

1.19 

0.996 

5.01 

7L-a 

16.2 

16.8 

2.64 

0.95 

1.29 

0.999 

5.46 

7L-b 

14.9 

22.5 

2.66 

1.19 

1.61 

0.998 

4.91 

7R-c 

17.7 

17.4 

1.77 

2.52 

2.07 

0.997 

5.83 

7R-d 

21.0 

8.3 

1.10 

0.85 

1.79 

0.999 

7.07 

7R-e 

24.7 

13.7 

1.37 

1.06 

1.62 

0.997 

8.69 

8R-a 

23.8 

26.6 

1.78 

2.58 

1.40 

0.992 

7.76 

8R-b 

22.2 

29.9 

1.44 

1.95 

1.11 

0.989 

6.09 

Mean 

17.5 

17.8 

1.67 

1.44 

1.50 

0.995 

5.80 

St.  Dev. 

4.2 

7.8 

0.71 

0.66 

0.30 

0.005 

1.33 
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TABLE  III.  Statistics  of  wavefront  distortion  caused  by  fourteen  tissue  paths  including  rib 
structures.  The  footnotes  associated  with  the  label  for  each  path  indicate  morphological  features 
and  physical  phenomena  that  affected  the  wavefront  distortion  computed  for  that  path.  The  format 
is  analogous  to  that  in  Table  II. 


Path 

Thickness 

(mm) 

ATF 

rms 

(ns) 

CL 

(mm) 

rms 

(dB) 

ELF 

CL 

(mm) 

WSF 

Attenuation 

(dB) 

4L-a1,2,3,4 

21.0 

260.3 

3.00 

2.58 

2.72 

0.968 

15.33 

4L-b2’3 

17.6 

161.9 

1.90 

4.16 

1.49 

0.641 

43.35 

5L-b2 

14.2 

92.5 

0.69 

3.06 

1.92 

0.775 

26.87 

7L-c3’5 

17.8 

47.2 

1.58 

5.33 

2.04 

0.958 

19.66 

7R-a1,3,4 

30.4 

123.1 

2.12 

3.80 

1.78 

0.960 

16.57 

7R-b3,5 

24.3 

165.6 

2.71 

6.88 

2.07 

0.274 

43.06 

8L-a3 

25.3 

113.9 

1.18 

7.75 

2.29 

0.907 

32.44 

SL-b1’4 

22.8 

109.7 

2.05 

3.43 

1.22 

0.974 

10.28 

8L-c3 

28.8 

134.0 

2.75 

3.04 

1.57 

0.944 

40.47 

8L-d4 

23.6 

78.9 

0.64 

3.06 

1.55 

0.950 

6.78 

8L-e3 

26.4 

208.8 

1.91 

3.62 

1.50 

0.810 

44.27 

8L-f1,4 

28.5 

169.9 

1.79 

5.02 

1.95 

0.916 

10.70 

8L-g3 

27.6 

210.8 

1.40 

3.36 

1.35 

0.892 

44.22 

8R-c2,3 

24.9 

81.4 

2.08 

2.76 

1.25 

0.962 

44.32 

1  Small  intercostal  spaces 
2Cancellous  bone 


3Cortical  bone  and  cartilage 

4Strong  rib  reflections 

5Cortical  bone  within  cartilage 
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FIGURES 


FIG.  1.  Chest  tissue  maps  used  in  simulations.  Simulated  apertures  are  indicated  using 
lower-case  letters  for  each  cross  section.  Smaller  arrows  indicate  55-element  (11.60  mm)  aper¬ 
tures  while  large  arrows  indicate  68-element  (14.28  mm)  apertures. 
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FIG.  2.  Summary  of  distortion  statistics  for  soft  tissue  paths.  The  bar  chart  shows  mean  values 
of  the  rms  arrival  time  fluctuations  (ATF),  rms  energy  level  fluctuations  (ELF),  correlations  lengths 
(CL)  of  these  fluctuations,  and  waveform  similarity  factors  (WSF)  for  the  simulations  performed 
in  the  present  paper  and  the  experiments  reported  in  Ref.  8.  Error  bars  indicate  a  range  of  plus 
or  minus  one  standard  deviation  from  the  mean. 
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V 


FIG.  3.  Measured  waveforms  for  three  propagation  paths  in  specimen  8L.  Each  panel  shows  re¬ 
ceived  waveforms  on  a  bipolar  logarithmic  gray  scale  with  a  dynamic  range  of  40  dB.  The  horizontal 
range  shown  in  each  panel  is  20  mm  and  the  vertical  range  shown  is  6.4  /zs.  (a)  Tissue  path  be¬ 
tween  two  ribs,  in  azimuth  direction  (parallel  to  ribs),  (b)  Path  including  a  rib,  azimuth  direction, 
(c)  Tissue  path  including  intercostal  space  between  two  ribs,  elevation  direction  (perpendicular  to 
ribs). 
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FIG.  4.  Propagation  through  the  central  rib  in  cross  section  5L  (path  5L-b).  Panels  (a)-(d) 
show  instantaneous  acoustic  pressure  fields  at  successive  intervals  of  2.17  /is.  Each  panel  shows  an 
area  that  spans  20.32  mm  horizontally  and  14.58  mm  vertically.  Logarithmically  compressed  wave- 
fronts  are  shown  on  a  bipolar  scale  with  black  representing  minimum  pressure,  white  representing 
maximum  pressure,  and  a  dynamic  range  of  57  dB. 
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FIG.  5.  Propagation  through  an  intercostal  space  in  cross  section  8L  (path  8L-b).  Panels 
(a)-(d)  show  instantaneous  wavefields  at  successive  intervals  of  3.62  ps.  Each  panel  shows  an  area 
that  spans  28.27  mm  horizontally  and  21.20  mm  vertically.  Wavefronts  are  shown  using  the  same 
format  as  in  Fig.  4. 


35 


30 


rms  ATF  (ns) 


rms  ELF  (dB) 


25  ■ 

20  -  t  _ 

15  - 
10  - 
5  - 

o  U — — U — . — U — — 

1.6  MHz  2.3  MHz  3.0  MHz 


CL-ATF  (mm)  CL-ELF  (mm) 


WSF 


FIG.  6.  Summary  of  frequency-dependent  distortion  results.  Mean  rms  arrival  time  fluctua¬ 
tions  (ATF),  energy  level  fluctuations  (ELF),  correlation  lengths  (CL)  of  these  fluctuations,  and 
waveform  similarity  factors  (WSF)  are  shown  for  each  of  the  three  pulse  frequencies  investigated. 
Error  bars  indicate  a  range  of  plus  or  minus  one  standard  deviation  from  the  mean. 
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Abstract 


Large-scale  simulation  of  ultrasonic  pulse  propagation  in  inhomogeneous  tissue  is  important 
for  study  of  ultrasound-tissue  interaction  as  well  as  for  development  of  new  imaging  methods. 
Typical  scales  of  interest  span  hundreds  of  wavelengths;  most  current  two-dimensional  methods, 
such  as  finite-difference  and  finite-element  methods,  are  unable  to  compute  propagation  on  this 
scale  with  the  efficiency  needed  for  imaging  studies.  Furthermore,  for  most  available  methods  of 
simulating  ultrasonic  propagation,  large-scale  three-dimensional  computations  of  ultrasonic  scat¬ 
tering  are  infeasible.  Some  of  these  difficulties  have  been  overcome  by  previous  pseudospectral 
and  A;-space  methods,  which  allow  substantial  portions  of  the  necessary  computations  to  be  exe¬ 
cuted  using  fast  Fourier  transforms.  This  paper  presents  a  new  A;-space  method  that  has  advantages 
of  both  past  fc-space  methods  and  pseudospectral  methods.  In  this  method,  the  spatial  differential 
equations  are  solved  by  a  simple  Fourier  transform  method  and  temporal  iteration  is  performed  us¬ 
ing  a  k-t  space  propagator.  The  applicability  of  the  A;-space  method  to  large-scale  soft-tissue  mod¬ 
eling  is  shown  by  simulating  propagation  through  several  tissue-mimicking  cylinders  as  well  as  a 
model  chest  wall  cross  section.  Numerical  results  indicate  that  this  method  is  accurate  for  large- 
scale  soft-tissue  computations,  with  much  greater  efficiency  than  that  of  an  analogous  leapfrog 
pseudospectral  method  or  a  2-4  finite  difference  time-domain  method.  However,  numerical  results 
also  indicate  that  the  /c-space  method  is  less  accurate  than  the  finite-difference  method  for  a  high- 
contrast  scatterer  with  bone-like  properties,  although  qualitatively  reasonable  results  can  still  be 
obtained  by  the  &-space  method  with  high  efficiency.  Possible  extensions  to  the  method,  including 
three-dimensional  computations  and  inclusion  of  absorption  effects,  are  discussed. 
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I.  INTRODUCTION 


Computation  of  a  scattered  acoustic  field,  given  an  incident  wavefield  and  complete  specification  of 
an  inhomogeneous  medium,  is  known  as  the  forward  scattering  problem.  In  contrast,  computation 
of  medium  parameters  given  the  incident  and  the  scattered  field  is  known  as  the  inverse  scattering 
problem.  Numerical  solution  of  the  forward  scattering  problem  is  central  to  many  aspects  of  ul¬ 
trasonic  imaging,  including  inverse  scattering  methods,  numerical  studies  of  wavefront  distortion, 
and  development  of  new  methods  for  adaptive  focusing. 

Most  methods  for  numerical  solution  of  the  forward  scattering  problem  fall  into  one  of  three 
categories:  finite-difference  methods,  finite-element  methods,  and  spectral  methods.  Finite- 
difference  and  finite-element  methods  are  known  as  local  because  the  wave  propagation  equations 
of  interest  are  solved  at  each  point  based  only  on  the  states  of  nearby  points.  In  contrast,  spectral 
methods  such  as  the  A; -space  method  [1]— [6]  and  the  pseudospectral  approach  [7]— [13]  are  called 
global  because  information  from  the  entire  wavefield  is  employed  to  solve  the  wave  propagation 
equations  at  each  point.  Because  of  their  global  nature,  spectral  methods  can  be  more  accurate 
than  local  methods — for  instance,  pseudospectral  methods  applied  to  periodic  problems  have  been 
shown  to  be  equivalent  to  finite-difference  methods  of  infinite  order  [11]. 

Spectral  methods  also  have  considerable  advantages  for  large-scale  forward  solvers  because  the 
required  storage  and  the  number  of  operations  per  iteration  can  be  dramatically  reduced  compared 
to  local  methods.  This  advantage  occurs  for  two  reasons.  First,  since  spectral  methods  contain 
derivatives  and/or  convolutions  that  can  be  evaluated  using  fast  Fourier  transforms,  the  number 
of  operations  per  iteration  is  proportional  to  N  log2  N  for  two-dimensional  computations,  while 
the  number  of  operations  per  iteration  is  proportional  to  N 2  for  implicit  finite-element  methods 
[14].  Second,  the  increased  accuracy  of  spectral  methods  can  allow  computations  to  be  performed 
on  coarser  grids  while  maintaining  accuracy.  For  example,  finite-element  methods  and  high-order 
finite-difference  methods  typically  require  grid  spacings  on  the  order  of  eight  points  per  minimum 
wavelength,  while  second-order  finite-difference  methods  can  require  twenty  points  per  wave¬ 
length  [9].  Spectral  methods,  in  theory,  require  only  two  points  per  wavelength  (spatial  Nyquist 
sampling),  although  for  computations  of  propagation  in  inhomogeneous  media,  greater  accuracy 
is  achieved  with  grid  spacings  on  the  order  of  four  points  per  wavelength  [9, 10, 13]. 

This  report  addresses  the  problem  of  large-scale  ultrasonic  wave  propagation  in  biological  me¬ 
dia  such  as  human  tissue.  For  problems  of  interest  in  medical  ultrasound,  domain  sizes  can  often 
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exceed  the  capabilities  of  conventional  forward  solvers.  For  example,  one  calculation  of  realistic 
scale  would  be  the  simulated  propagation  of  a  pulse  with  an  upper  bandwidth  limit  of  5  MHz  in 
a  volume  of  dimensions  30  mm  on  each  side  and  a  nominal  sound  speed  of  1.5  mm///s,  so  that 
the  minimum  wavelength  is  0.3  mm.  For  this  calculation,  a  second-order  finite-difference  method 
(using  twenty  points  per  wavelength)  would  require  a  three-dimensional  grid  containing  8  x  1010 
nodes,  a  finite-element  or  fourth-order  finite-difference  method  (using  ten  points  per  wavelength) 
would  require  1  x  109  nodes,  and  a  spectral  method  (using  four  points  per  wavelength)  would  re¬ 
quire  6.4  x  107  nodes.  Since  a  grid  of  6.4  x  107  single-precision  complex  numbers  requires  storage 
of  512  megabytes,  only  spectral  methods  are  feasible  for  realistic  three-dimensional  computa¬ 
tions  on  present-day  computers  that  typically  have  a  maximum  random-access  memory  storage  of 
several  gigabytes.  The  efficiency  provided  by  fast  Fourier  transform  implementations  of  spectral 
algorithms  is  a  further  reason  why  spectral  methods  are  a  practical  approach  to  large-scale  and 
three-dimensional  computations  of  ultrasonic  wave  propagation. 

Previous  spectral  approaches  have  included  pseudospectral  methods,  in  which  spatial  deriva¬ 
tives  are  evaluated  globally  using  Fourier  transformation  and  wavefields  are  advanced  in  time  using 
various  numerical  integration  techniques  [7]— [13].  This  method  has  provided  high  accuracy  in  may 
cases;  however,  temporal  iteration  techniques  that  provide  good  accuracy  for  large-scale  models 
typically  require  significant  additional  computations  and/or  storage  of  wave  fields  from  additional 
time  steps  [12],  so  that  the  efficiency  advantages  of  the  pseudospectral  approach  are  less  than  might 
first  be  expected.  The  fc-space  family  of  methods  [1]— [6]  can  overcome  this  problem  by  providing 
explicit  temporal  propagators  related  to  the  Green’s  function  for  wave  propagation  in  k-t  (spatial 
frequency  and  time)  space.  Previous  implementations  of  the  fc-space  method  have  posed  the  spa¬ 
tial  part  of  the  wave  propagation  equations  in  an  integral  form  that  can  be  solved  using  discrete 
convolution  [l]-[6].  Although  this  method  has  been  shown  to  be  efficient  and  accurate,  the  dis¬ 
crete  convolution  step,  which  employs  inherently  singular  spatial  Green’s  functions,  can  present 
complexities  in  numerical  implementations. 

The  present  paper  presents  an  alternative  derivation  of  the  A;-space  method  using  a  differential 
representation  of  the  wave  propagation  equations.  The  spatial  part  of  the  wave  propagation  equa¬ 
tions  is  solved  using  Fourier  transforms  in  a  manner  analogous  to  past  pseudospectral  methods; 
this  derivation  is  shown  to  be  theoretically  equivalent  to  previous  integral  formulations  of  the  k- 
space  method.  Temporal  iteration  is  performed  using  a  k-t  space  propagator  [2],  which  is  shown  to 
provide  much  greater  accuracy  than  “leapfrog”  iteration  without  significant  additional  computation 
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or  storage  requirements.  Thus,  the  new  fc-space  method  has  advantages  previously  associated  with 
both  pseudospectral  and  integral  k- space  methods. 

Below,  a  derivation  of  the  new  A;-space  method  is  presented  for  propagation  in  fluid  media 
with  spatially-dependent  sound  speed  and  density.  For  several  canonical  forward  problems  rel¬ 
evant  to  ultrasonic  imaging,  the  accuracy  and  efficiency  of  the  £;-space  method  is  compared  to 
a  pseudospectral  method  employing  leapfrog  iteration  and  to  a  2-4  finite  difference  time-domain 
method.  The  k- space  and  finite-difference  methods  are  also  used  in  an  example  computation  for 
a  large-scale  two-dimensional  tissue  model.  Possible  extensions  of  the  new  A;-space  method,  in¬ 
cluding  multiple  relaxation  effects  for  absorption,  full  elastic  wave  propagation,  perfectly  matched 
layers  for  absorbing  boundary  conditions,  three-dimensional  computations,  and  parallelization,  are 
discussed. 
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II.  THEORY 

A.  Derivation  of  the  fc-space  method 


A  new  k-space  method  for  solving  the  acoustic  forward  scattering  problem  is  derived  here.  The 
method  is  applicable  to  large-scale  modeling  of  linear  ultrasonic  propagation  in  soft  tissues,  which 
are  modeled  here  as  fluid  media  with  spatially-dependent  sound  speed  and  density.  Although  the  k- 
space  method  described  below  can  be  extended  to  include  absorption  effects,  acoustic  nonlinearity, 
and  shear-wave  propagation,  these  effects  are  neglected  in  this  derivation  for  simplicity. 

For  a  fluid  medium  with  spatially-dependent  sound  speed  and  density,  the  wave  equation  is 


1 


,P(X) 


Vp(x,  t) 


d2p(x,  t) 


0, 


(1) 


p(x)  c(x)2  dt2 

where  p(x,  t)  is  the  acoustic  perturbation  in  pressure,  p(x)  is  the  ambient  density,  and  c(x)  is  the 
ambient  sound  speed. 

By  defining  the  normalized  wavefield 


/(x,f)  =p(x,t)/v/p(x,t), 


Eq.  (1)  can  be  rewritten  in  the  form 

tin-*  1  aV(x,i) 


V2/(x,«) 


c2(x)  dt2 


-  K*)f  Cm)  =  o, 


(2) 


(3) 


where 


&(x)  =  \fp (x)  v2  (l/ yp(x)^ .  (4) 

The  normalized  field  /  may  now  be  split  into  an  incident  and  a  scattered  part,  /(x,  t)  = 
fi(x,  t)  +  /s(x,  t).  The  normalized  incident  field  ft  satisfies  the  homogeneous  wave  equation 


*/,<*,«) o, 


(5) 


Cq  dt 2 

where  c0  is  a  nominal  “background”  sound  speed  for  the  inhomogeneous  medium,  while  the  nor¬ 
malized  scattered  field  fs  satisfies  the  inhomogeneous  equation 


V2fs{x,t)  - 


1  92/a(x,f) 

c2  dt2 


,c(x)5 


co, 


<92/(x,*) 

dt 2 


+  b(x)f{x,t). 


(6) 


The  terms  on  the  right  hand  side  of  Eq.  (6)  are  effective  sources  of  scattered  waves,  and  can  be 
written  in  terms  of  a  source  associated  with  sound  speed  variations, 


v(x,t)  = 


r2 

co 


-i  /(x>*)> 


(7) 
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and  a  source  associated  with  density  variations. 


=  c20b(x)  f(x,t).  (8) 

With  the  additional  definition  of  an  auxiliary  field  w(x,t)  =  fs(x,  t)  +  v(x,t),  Eq.  (6)  can  be 
written  as 

V2/,(x, t)  =  -  +  6(x) /(x, t).  (9) 

After  spatial  Fourier  transformation  of  Eq.  (9),  the  inhomogeneous-medium  wave  equation  can  be 
written  as  the  coupled  set  of  equations 

d2^(2M)  =  (c0k)2  [V(k,  t)  -  W{ k,  t)]  -  Q( k,  t)  (10) 

v(x,t)  =  d(x)  [/i(x,f)  +  tu(x,i)]  (11) 

q(x,t)  =  c2  b(x)  [fi(x,  t)  +  w(x,  t)  —  v(x,  t)],  (12) 

where  capital  letters  denote  spatially  Fourier  transformed  variables  and  the  sound  speed  contrast 

function  d(x)  is  defined  as 

d(x)  =  1  -  «W!.  (13) 

co 

In  numerical  implementation  of  the  k- space  algorithm,  Eq.  (10)  is  used  to  advance  the  auxiliary 
field  J^(k,  t)  in  time,  Eq.  (11)  updates  the  effective  source  v  associated  with  sound  speed  varia¬ 
tions,  and  Eq.  (12)  updates  the  effective  source  q  associated  with  density  inhomogeneities.  Notable 
is  that  the  effective  source  v  is  directly  proportional  to  the  square  of  the  sound  speed  variation  of 
the  medium,  while  the  effective  source  q  is  directly  proportional  to  the  Laplacian  of  \/^J p{x). 
Thus,  for  a  piecewise-constant  inhomogeneous  medium,  v  may  be  non-zero  everywhere  while  q  is 
nonzero  (and  singular)  only  on  borders  between  regions. 

The  algorithm  derived  above  can  be  considered  an  extension  of  Bojarski’s  formulation  [1,  2], 
which  was  obtained  using  an  integral  representation  of  the  acoustic  wave  equation.  This  can  be 
shown  by  rewriting  Eq.  (10)  using  the  k-t  space  Green’s  function  [2] 

G(k,  t-r)  =  c0k  sin[coA;(t  -  t)]  H(t  -  r),  (14) 

where  H  is  the  Heaviside  step  function.  The  Green’s  function  defined  by  Eq.  (14)  satisfies  the 
differential  equation 

d2G^t]  +  c2k2G(k,t)  =  c2k26(t).  (15) 
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dr. 


(16) 


For  zero  initial  conditions,  Eq.  (10)  can  then  be  rewritten  as 

W (k,  t)  =  c0k  sin[c0A:(t  -  r)]  (k,  r)  -  k,  r) 

For  a  medium  of  constant  density,  Q  =  0  and  Eq.  (16)  reduces  to  Eq.  (56)  of  Ref.  [2],  Thus,  the 
present  method  can  be  considered  a  generalized  form  of  Bojarski’s  k- space  method. 

Equations  (10)— (12)  can  be  written  in  several  different  discrete  forms.  For  instance,  using  a 
second-order  accurate  finite-difference  representation  of  the  second-order  time  derivative,  Eq.  (10) 
can  be  written  in  the  discrete  form 

V(k,t)-W(k,t)-^0  .  (17) 

where  A t  is  the  time  step.  This  is  known  as  “leapfrog”  iteration;  use  of  Eq.  (17)  in  the  cur¬ 
rent  method  is  analogous  to  commonly  used  pseudospectral  approaches  [12,  13].  (Although  in¬ 
creased  accuracy  can  be  achieved  by  higher-order  methods  such  as  fourth-order  Adams-Bashforth 
or  Adams-Moulton  iteration,  these  methods  have  the  disadvantage  of  requiring  storage  of  the 
entire  computational  grid  for  additional  time  steps  [11,  12].)  A  more  accurate  form,  how¬ 
ever,  results  from  expanding  Eq.  (16)  about  the  instants  t  ±  At  and  making  the  approximation 
V(k,t)  «  (l/2)[V(k,t  +  At/2)  +  V{k,t  -  At/2 )]  [2].  This  yields  the  k-t  space  propagation 
equation 

W(k,t  +  At)-2W(k,t)  +  W(k,t-At)=4sm2h^]  V(k,t)  -  W(k,t)  -  ^0  . 

(18) 

Alternatively,  Eq.  (18)  can  be  derived  from  Eq.  (10)  using  the  method  of  mode-dependent 
discretization  [15],  in  which  a  difference  operator  is  chosen  to  produce  the  same  result  as  the  dif¬ 
ferential  operator  on  a  set  of  chosen  “coincident”  functions.  A  natural  choice  for  these  coincident 
functions  is  the  set  of  null  functions  for  the  differential  operator  of  interest.  Here,  these  are  e±lCokt 
for  the  operator  d2/dt 2  -I-  ( c0k )2  that  is  applied  to  W  in  Eqs.  (10)— (12).  By  applying  the  differential 
operator  d2/dt 2  +  ( c0k )2  as  well  as  its  discrete  form  to  these  coincident  functions  and  equalizing 
the  results,  Eq.  (18)  follows.  It  may  also  be  noted  that  Eqs.  (17)  and  (18)  are  equivalent  in  the  limit 
of  small  At.  However,  results  shown  below  indicate  that  use  of  the  k-t  propagator  (18)  provides 
much  greater  accuracy  for  larger  time  steps. 

The  present  /c-space  algorithm  can  now  be  summarized  as  follows: 


W(k,  t  +  At)-2  IF(k,  t)  +  W( k,  t- At)  =  (c0kAt)2 


1.  Define  the  contrast  functions  d(x)  and  6(x). 


2.  Initialize  the  field  w(x,  t )  to  zero. 

3.  Define  the  incident  wave  /4(x,  t)  on  the  entire  grid. 

4.  Compute  v(x,  t)  using  Eq.  (11). 

5.  Compute  q(x,  t )  using  Eq.  (12). 

6.  Obtain  V"(k,  t),  Q(k,  t),  and  VF(k,f)  by  FFT. 

7.  Evaluate  W (k,  t  +  A t)  from  Eq.  (17)  or  Eq.  (18). 

8.  Obtain  w(x,  t  +  At)  by  inverse  FFT. 

9.  Set  t  -¥  t  +  At  and  go  to  (3). 

To  distinguish  between  the  standard  leapfrog  iteration  method  and  the  improved  method  used 
here,  the  following  nomenclature  is  used  in  the  present  paper.  The  above  algorithm  employing 
Eq.  (17)  for  temporal  iteration  is  referred  to  as  a  leapfrog  pseudospectral  method  below,  while 
the  algorithm  employing  Eq.  (18)  is  referred  to  as  a  /c-space  method.  This  nomenclature  is  used 
because,  as  shown  above,  the  algorithm  employing  Eq.  (18)  is  equivalent  to  an  extended  form  of 
Bojarski’s  k- space  method  [2]  cast  in  terms  of  differential  equations  rather  than  integral  equations. 

B.  Numerical  stability  and  accuracy 

The  stability  of  the  method  derived  above  can  be  evaluated  using  standard,  linear  von  Neumann 
stability  analysis  [16].  In  this  technique,  the  difference  equations  that  comprise  Eqs.  (17)  and  (18) 
are  applied  to  a  test  function 

Wiest(k,nAt)  =  (f>(k)n  tp(k),  (19) 

where  k)  is  a  spatial-frequency  domain  eigenmode  and  <j>( k)  is  a  temporal  amplification  factor. 
If  a  difference  equation  admits  solutions  with  |0(k)  |  >  1  for  any  vector  wavenumber  k,  errors  may 
grow  exponentially  with  time  and  the  solution  is  thus  unstable.  If  |0(k)  |  <  1  for  all  wavenumbers, 
then  the  solution  is  numerically  stable.  For  simplicity,  since  results  can  depend  on  the  inhomogene¬ 
ity  employed,  this  stability  computation  is  performed  here  under  the  assumption  of  a  homogeneous 
medium. 
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Application  of  this  technique  to  Eq.  (17),  which  represents  a  leapfrog  pseudospectral  approach, 
yields  a  quadratic  equation  for  (j>( k).  The  resulting  stability  condition  is 

Cmax^maxAf  <  2,  (20) 


where  cmax  is  the  maximum  sound  speed  in  the  region  of  computation,  kmax  =  n/Ax  is  the  max¬ 
imum  wavenumber  in  the  discrete  Fourier  transforms  used  to  compute  W (k,  t),  and  At  and  Ax, 
respectively,  are  the  temporal  and  spatial  steps  employed.  Using  the  standard  definition  for  a 
Courant-Friedrichs-Lewy  (CFL)  number  [17],  the  stability  condition 

CFL  =  ^  <  --25-  (21) 

/\X  7T  Cmax 


is  obtained  for  the  leapfrog  pseudospectral  method  represented  by  Eq.  (17). 

Application  of  the  same  analysis  to  the  /c-space  method  represented  by  Eq.  (18)  yields  a  very 
different  result:  in  the  absence  of  inhomogeneities,  the  linear  numerical  stability  of  the  fc-space 
method  is  unconditional.  An  upper  limit  on  the  minimum  required  time  step,  however,  arises  from 
the  requirement  of  sampling  at  the  Nyquist  rate:  that  is,  the  time  step  should  be  sufficiently  small 
to  allow  two  samples  per  period  for  the  highest-frequency  component  of  the  computed  field.  Thus, 
the  requirement  can  be  written 


At  < 


7 T 


Ax 


2/n 


-max  ^  max 


(22) 


or  simply  CFL  <  c0/cmax. 

The  Fourier  transform  steps  in  the  above  algorithm  can  lead  to  numerical  artifacts  (related  to  the 
Gibbs  phenomenon)  when  the  scattering  object  contains  discontinuities  in  sound  speed  or  density. 
To  avoid  such  artifacts,  the  scattering  object  can  be  spatially  filtered  to  smooth  discontinuities. 
That  is,  the  sound-speed  contrast  function  c20{x)/c2  - 1  and  the  density  contrast  function  of  Eq.  (4) 
are  replaced  by  filtered  functions  of  the  form 


^filtered  (x)  =  F  —  1[[/(k)  </>(k)], 


(23) 


in  which  the  Fourier  transform  U{ k)  of  the  function  u(x)  is  multiplied  by  a  low-pass  spatial- 
frequency  filter  (f>( k).  In  the  present  study,  the  filter  employed  is  the  spatial-frequency  Blackman 
window  [18] 

(7T fc  \  /  2ttAj  \ 

- — )  +  0.08  cos  ( - — ).  (24) 

max  /  '  ^max  / 
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Spatial-frequency  windowing,  like  that  represented  by  Eq.  (24),  naturally  reduces  high  spatial- 
frequency  components  of  the  wavefield,  so  such  windowing  may  also  allow  use  of  larger  time  steps 
while  still  meeting  the  Nyquist  criterion  of  Eq.  (22). 
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III.  Numerical  Methods 


Numerical  implementation  of  the  new  k- space  algorithm  was  accomplished  using  the  algorithm 
described  above.  The  normalized  incident  wave  /*(x,  t)  was  defined  as  a  plane  wave  with  Gaussian 
temporal  shape, 

/j(x,  t)  =  p(x)“5  sin(o;oT)  e-r2/(2£r2),  (25) 

where  r  is  the  retarded  time  r  =  t  -  (x  -  x0)/c0  and  x0  is  the  initial  central  position  of  the  wave. 
Boundary  conditions  were  implicitly  periodic. 

Wavefields  were  computed  on  two-dimensional  grids  chosen  to  be  large  enough  to  avoid  in¬ 
fluence  of  “wraparound”  error  within  the  temporal  window  of  interest.  All  k- space  computations 
were  performed  on  square  grids  of  size  N  by  N.  Prior  to  execution  of  the  main  computation  loop, 
the  Laplacian  occurring  in  Eq.  (4)  was  evaluated  using  second-order  accurate,  centered  finite- 
difference  representations  of  the  second  derivative  in  each  direction.  Within  the  main  computa¬ 
tional  loop,  all  spatial  derivatives  were  evaluated  by  Fourier  transformation,  implemented  using  a 
fast  Fourier  transform  (FFT)  algorithm  [20].  For  maximum  FFT  efficiency,  grid  sizes  N  were  cho¬ 
sen  to  be  integers  with  prime  factors  no  higher  than  3.  In  some  cases,  wavefields  were  computed 
on  a  grid  of  size  smaller  than  N  and  the  fields  were  zero-padded  to  size  N  by  N. 

The  spatial-frequency  time-domain  wavefield  W  (k,  t  +  At)  was  windowed  using  the  rectan¬ 
gular  window 

<K k)  =  H(kmax  -  k),  (26) 

before  inversion  to  yield  w{x,  t  +  At).  (That  is,  between  steps  7  and  8  in  the  algorithm  enumerated 
above.)  In  Eq.  (26),  H  is,  as  before,  the  Heaviside  step  function,  kmax  is  the  maximum  wavenumber 
magnitude  (equal  to  'k/ Ax),  and  k  is  the  magnitude  of  the  vector  wavenumber  k.  In  some  cases, 
the  contrast  functions  6(x)  and  d(x)  were  also  smoothed  by  windowing  in  the  spatial-frequency 
domain  using  Eq.  (24)  with  a  wavenumber  cutoff  equal  to  7r/ Ax. 

For  comparison,  wavefields  were  also  computed  using  a  second-order  in  time,  fourth-order  in 
space  finite-difference  method,  described  in  Refs.  [17]  and  [19]-[23].  For  the  finite-difference 
computations,  the  incident  wave  was  specified  by  a  single  initial  condition  rather  than  updated 
at  each  time  step.  Periodic  boundary  conditions  were  applied  on  the  sides  perpendicular  to  the 
wavefront,  while  first-order  radiation  boundary  conditions  [23]  were  applied  on  the  sides  parallel  to 
the  wavefront.  Time  steps  were  determined  using  a  CFL  number  of  0.25,  which  is  a  natural  choice 
for  this  finite-difference  method  [21].  To  reduce  the  computational  burden,  domain  sizes  were 
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minimized  in  the  dimension  for  which  the  radiation  boundary  conditions  eliminated  wraparound 
errors.  As  in  Refs.  [23]— [25],  computations  were  performed  at  each  time  step  only  on  portions  of 
the  grid  where  the  wavefields  were  nonzero;  this  reduces  the  required  computation  time  for  the 
finite-difference  method  by  about  one  half. 

To  quantitatively  test  the  accuracy  of  the  A; -space  and  finite-difference  methods,  benchmark 
computations  were  performed  using  an  exact  series  solution  for  the  scattering  of  a  plane  wave  by  a 
fluid  cylinder  [26].  The  sampling  rate  and  waveform  shape  were  chosen  to  match  the  time-domain 
simulation  data  for  the  case  of  interest.  The  pressure  field  was  then  computed  for  each  frequency 
component  with  relative  magnitude  within  60  dB  of  the  peak  magnitude.  Each  single-frequency 
computation  truncated  the  series  at  the  term  having  a  relative  contribution  less  than  10~12  times  the 
sum  of  all  terms.  The  frequency-domain  scattered  fields  were  then  inverted  by  FFT  to  obtain  exact 
numerical  solutions  for  the  time-domain  pressure  fields  at  the  simulated  measurement  points. 

Benchmark  studies  of  accuracy  were  performed  using  a  cylinder  with  radius  2.0  mm  and  acous¬ 
tic  properties  of  human  fat,  and  a  background  medium  with  acoustic  properties  of  water  at  body 
temperature.  Rationale  for  use  of  these  values  is  discussed  in  Ref.  [23].  The  cylinder  had  a  sound 
speed  of  1.478  mm/^s  and  a  density  of  0.950  g/mm3,  while  the  background  medium  had  a  sound 
speed  of  1.524  mm//is  and  a  density  of  0.993  g/mm3.  The  incident  pulse  was  a  plane  wave  with 
Gaussian  temporal  characteristics,  a  temporal  Gaussian  parameter  a  =  0.25  //s,  and  a  central 
starting  position  of  x  =  -4.5  mm  at  time  zero.  For  this  pulse,  a  nominal  maximum  frequency 
is  4.43  MHz,  corresponding  to  the  spectral  point  40  dB  down  from  the  center  frequency  (for  the 
benchmark  problem,  this  frequency  corresponds  to  a  minimum  wavelength  of  0.334  mm).  The 
fc-space,  leapfrog  pseudospectral,  finite-difference,  and  exact  methods  described  above  were  used 
to  compute  time  histories  of  the  total  pressure  field  at  128  equally-spaced  “measurement”  points 
spanning  a  circle  of  radius  2.5  mm  concentric  to  the  cylinder.  The  domain  size  for  each  computa¬ 
tion  employing  this  cylinder  was  9.5  x  9.5  mm2. 

Further  studies  of  accuracy  were  performed  using  a  cylinder  of  radius  10  mm.  Other  parameters 
were  as  described  above  for  the  small  problem,  except  that  the  radius  of  the  measurement  circle  was 
12.5  mm  and  the  starting  position  of  the  wavefront  was  x  =  —14.5  mm.  The  k- space  method  was 
employed  to  compute  two  cases  corresponding  to  unsmoothed  and  smoothed  contrast  functions, 
using  a  spatial  step  of  four  points  per  minimum  wavelength  and  a  CFL  number  of  0.5.  In  each 
A;-space  computation  for  this  cylinder,  the  domain  size  employed  was  72  x  72  mm2.  The  finite- 
difference  method  was  employed  to  compute  a  single  case,  using  a  spatial  step  of  fourteen  points 
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per  minimum  wavelength  and  a  CEL  number  of  0.25.  Because  storage  requirements  proved  to 
be  large  for  this  finite-difference  computation,  the  computational  configuration  was  modified  to 
employ  a  finite-width  plane  wave  with  a  13%  cosine  taper,  as  well  as  radiation  conditions  on  all 
computational  boundaries.  These  modifications,  which  were  performed  only  for  the  large  cylinder 
case,  allowed  reduction  of  the  domain  size  to  32  x  32  mm2,  so  that  the  computation  could  be 
performed  within  the  storage  limits  of  the  computer  employed. 

To  evaluate  the  relative  accuracy  and  efficiency  of  the  fc-space  and  finite-difference  methods 
for  a  high-contrast  scatterer,  computations  were  also  performed  using  a  cylinder  of  radius  2.0  mm 
with  the  sound  speed  and  density  of  human  bone.  The  values  employed  were  a  sound  speed 
of  3.54  mm/yus  and  a  density  of  1.99  g/mm3,  as  in  Ref.  [25].  The  incident  pulse,  receiver,  and 
computational  domain  characteristics  were  identical  to  those  for  the  2.0  mm  “fat”  cylinder  case 
described  above. 

In  all  of  the  above  accuracy  tests,  a  quantitative  measure  of  the  accuracy  was  obtained  using 
the  time-domain  L2  error  of  each  approximate  pressure  field  papProx(x,  t)  versus  the  corresponding 
exact  series  solution  pexact(x,  t).  This  quantity  has  the  definition 

_  ff  [^approx  (xr  1 1)  Pex act  fop  t)|  dxr  dt  (27) 

If  IPexact (Xr ,t)\2dxrdt 

where  the  integrals  are  evaluated  using  discrete  summation  over  all  receiver  points  xr  and  the  en¬ 
tire  time  span  of  the  computations.  Eq.  (27)  represents  an  accuracy  criterion  that  is  much  stricter 
than  more  general  criteria,  such  as  comparison  of  the  rms  waveform  amplitude  or  the  amplitude 
and  phase  at  the  center  frequency.  To  achieve  a  low  L2  error  by  the  definition  of  Eq.  (27),  both  the 
waveform  amplitude  and  phase  must  be  accurately  computed  for  all  significant  frequency  compo¬ 
nents  of  the  field. 

Lastly,  the  use  of  the  present  fc-space  method  in  a  more  realistic  simulation  of  ultrasonic  propa¬ 
gation  was  also  tested.  For  this  purpose,  a  cross-sectional  tissue  map  of  the  human  chest  wall  [25] 
was  used  as  the  simulated  medium.  A  pulse  center  frequency  of  3.0  MHz  was  employed  together 
with  a  temporal  Gaussian  parameter  of  0.4766  //s;  these  parameters  correspond  to  the  highest  cen¬ 
ter  frequency  employed  in  the  simulation  study  reported  in  Ref.  [25].  The  corresponding  nominal 
minimum  wavelength  is  0.3252  mm.  The  k- space  computation  employed  4  points  per  minimum 
wavelength,  a  CFL  number  of  0.5,  and  a  grid  size  of  54.9  x  54.9  mm2.  The  finite-difference  com¬ 
putation  employed  14  points  per  minimum  wavelength,  a  CFL  number  of  0.25,  and  a  grid  size  of 
38.5x29.7  mm2. 
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IV.  Numerical  Results 


An  example  /c-space  computation,  performed  using  the  2.0  mm  cylinder  with  acoustic  properties 
of  human  fat,  is  illustrated  in  Fig.  1.  The  cylinder  boundary  is  also  sketched  in  each  panel.  For 
the  computation  shown,  the  contrast  functions  d(x)  and  6(x),  which  are  related  to  the  spatially- 
dependent  sound  speed  and  density  by  Eqs.  (13)  and  (4),  respectively,  were  smoothed  using  Eq.  24. 
The  time  history  of  the  total  wavefield  is  shown  as  computed  by  the  A; -space  method  for  a  spatial 
step  size  of  four  points  per  minimum  wavelength  and  a  CFL  number  of  0.5.  Details  visible  include 
a  scattered  wave  from  the  edge  nearest  the  initial  wavefront  (panel  3),  some  spurious  oscillations 
near  the  -60  dB  level  (panels  4-6),  weak  focusing  near  the  trailing  edge  of  the  cylinder  (panel  5), 
and  scattering  from  the  trailing  edge  (panels  6-9). 

Results  of  accuracy  benchmarks  for  the  /c-space  and  leapfrog  pseudospectral  methods  described 
above  are  shown  in  Fig.  2.  Each  of  these  computations  was  made  using  the  2.0  mm  radius  cylinder 
described  above  and  a  spatial  step  size  of  four  points  per  maximum  wavelength.  The  results  show 
that  the  /c-space  method  employing  the  k-t  space  propagator  of  Eq.  (18)  provides  much  higher 
accuracy  than  the  pseudospectral  method  employing  the  leapfrog  propagator  of  Eq.  (17).  The  two 
methods  provide  equivalent  results  for  very  small  time  steps  (CFL  numbers  less  than  about  0.1),  but 
the  /c-space  method  maintains  its  highest  accuracy  up  to  a  CFL  number  of  about  0.7.  In  contrast, 
the  pseudospectral  method  rapidly  increases  in  error  for  CFL  numbers  above  0.1. 

Error  results  for  the  pseudospectral  computations  shown  in  Fig.  2  are  not  given  for  CFL  num¬ 
bers  above  0.6  because,  for  higher  CFL  numbers,  the  computation  was  unstable.  (Computed  fields 
incurred  spurious  exponential  growth,  resulting  in  numerical  overflow.)  This  observation  of  in¬ 
stability  is  consistent  with  the  theoretical  stability  limit  of  0.6366  given  by  Eq.  (21)  for  this  case. 
Although  the  /c-space  method  did  not  incur  any  numerical  instability  for  the  range  of  CFL  num¬ 
bers  investigated,  the  error  of  this  method  grows  as  the  CFL  number  approaches  and  exceeds  unity, 
consistent  with  the  Nyquist  sampling  criterion  given  by  Eq.  (22). 

The  relative  accuracy  of  the  /c-space  and  finite-difference  methods  are  compared  in  Fig.  3  as  a 
function  of  the  spatial  step  size.  For  these  computations,  the  CFL  number  of  the  /c-space  compu¬ 
tations  was  held  constant  at  0.5,  consistent  with  the  range  of  accuracy  shown  in  Fig.  2,  while  the 
CFL  number  of  the  finite-difference  computations  was  held  at  0.25  [21],  Both  methods  achieve 
high  accuracy  as  the  grid  spacing  becomes  finer;  however,  the  k- space  method  achieves  higher  ac¬ 
curacy  for  much  larger  spatial  step  sizes.  The  L2  error  drops  below  0.05  for  k- space  computations 


14 


employing  only  four  points  per  minimum  wavelength,  while  achievement  of  the  same  accuracy 
criterion  requires  16  points  per  minimum  wavelength  for  the  finite-difference  computations.  This 
difference  suggests  that  storage  requirements  for  fc-space  computations  can  be  much  smaller  than 
those  for  finite-difference  computations  of  comparable  accuracy:  on  the  order  of  16  times  smaller 
for  two-dimensional  computations  and  64  times  smaller  for  three-dimensional  computations. 

Visual  comparison  of  simulated  waveforms  for  the  2.0  mm  diameter  cylinder  is  shown  in  Fig.  4. 
Waveforms  in  this  figure  are  those  computed  using  the  fc-space  (four  points  per  minimum  wave¬ 
length,  CFL  number  0.5,  with  both  unsmoothed  and  smoothed  contrast  functions),  finite-difference 
time-domain  (14  points  per  minimum  wavelength,  CFL  number  0.25),  and  exact  methods.  The  k- 
space  solution  for  the  unsmoothed  cylinder  shows  the  smallest  time-domain  L2  error  (0.0347),  but 
also  exhibits  spurious  waves  (nearly  60  dB  down  from  the  peak  pressure  amplitude)  between  the 
two  main  arrivals.  These  spurious  waves  are  removed  by  use  of  the  &-space  method  with  smoothed 
contrast  functions  (i.e.,  6(x)  from  Eq.  (4)  and  d(x)  from  Eq.  (13)  smoothed  using  Eq.  (24)  with 
kmax  —  7t/Ax),  but  the  L2  error  is  increased  to  0.0414  by  this  smoothing.  Visual  inspection  of 
Fig.  4(b)  indicates  that  the  amplitude  of  the  backscattered  wave  is  underpredicted  by  the  /c -space 
method  when  smoothed  contrast  functions  are  employed.  The  finite-difference  result  bears  a  strong 
qualitative  resemblance  to  the  exact  solution,  but  the  larger  L2  error  (0.0529)  indicates  that  phase 
errors  have  been  introduced  by  the  dispersion  inherent  to  the  finite-difference  method.  Computa¬ 
tion  times  (CPU  times  for  a  Linux  workstation  with  an  AMD  K6  processor  running  at  200  MHz) 
were  4.15  minutes  for  the  k- space  method  and  46.8  minutes  for  the  finite-difference  method,  so 
that  the  fc-space  method  yields  comparable  accuracy  at  much  less  computational  cost. 

Waveforms  for  the  10  mm  diameter  cylinder  are  shown  in  Fig.  5  in  a  format  analogous  to 
that  of  Fig.  4.  These  results  indicate  that,  as  for  the  smaller  cylinder,  smoothing  of  the  contrast 
functions  produces  a  reduction  in  spurious  low-amplitude  waves.  For  this  problem,  unlike  the 
2.0  mm  diameter  cylinder  discussed  above,  this  smoothing  slightly  increases  the  overall  accuracy. 
(The  time-domain  L2  error  is  0.1080  for  the  smoothed  case  versus  0.1103  for  the  unsmoothed 
case.)  The  finite-difference  solution,  using  14  points  per  wavelength  and  a  CFL  number  of  0.25, 
requires  longer  computation  time  (14.50  hours  vs.  4.61  hours)  and  produces  waveforms  with  poorer 
accuracy  (an  L2  error  of  0.1763)  than  the  fc-space  method. 

Results  for  the  2  mm  “bone”  cylinder  are  shown  in  Fig.  6.  In  this  case,  the  /c -space  method  using 
a  CFL  number  of  0.5  exhibited  numerical  instability,  possibly  associated  with  the  failure  of  this 
CFL  number  to  meet  the  Nyquist  condition  of  Eq.  (22).  To  maintain  a  consistent  temporal  sampling 
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rate,  the  time  step  was  reduced  in  proportion  to  the  increase  in  cmax,  resulting  in  a  CFL  number  of 
0.2153.  Required  computation  time  for  the  k- space  method  was  9.9  minutes;  the  time-domain  L2 
error  was  0.3072  for  the  unsmoothed  case  and  0.8551  for  the  smoothed  case.  As  can  be  seen  in 
Fig.  6,  the  k- space  computation  with  smoothing  reduces  the  level  of  spurious  wave-like  artifacts, 
but  shows  some  wavefield  details  less  accurately  than  the  computation  without  smoothing.  For 
the  k- space  method  without  smoothing,  reduction  of  the  time  step  by  an  additional  factor  of  two 
resulted  in  negligible  accuracy  improvement  (L2  error  0.3062)  and  a  doubling  of  the  computation 
time  to  19.7  minutes. 

For  the  “bone”  cylinder,  the  finite-difference  method,  employing  14  points  per  wavelength  and 
a  CFL  number  of  0.1076  (also  changed  in  proportion  to  cmax),  achieved  an  L2  error  of  0.0541  in  a 
computation  time  of  189  minutes.  This  result  indicates  that  finite-difference  methods  can  be  much 
more  accurate  than  k- space  methods  for  scattering  problems  involving  very  high-contrast  inhomo¬ 
geneities  such  as  bone  within  soft  tissue.  This  discrepancy  is  likely  to  be  associated  with  inaccurate 
representation  of  high-spatial  frequency  components  of  the  effective  sources  However,  the  fc-space 
solution,  as  seen  in  Fig.  6,  still  shows  good  qualitative  agreement  with  the  exact  solution. 

The  relative  inaccuracy  of  the  k- space  method  for  high-contrast  scatterers  may  be  associated 
with  aliasing  effects,  as  suggested  in  Ref.  [5].  That  is,  large  jumps  in  spatial  contrast  functions  are 
associated  with  significant  high-frequency  components  of  the  corresponding  A; -space  spectra.  If  the 
spatial-frequency  range  employed  in  the  A;-space  algorithm  is  not  sufficiently  large,  aliasing  errors 
result.  Smoothing  the  contrast  functions  removes  this  aliasing  but  introduces  additional  errors  be¬ 
cause  the  contrast  functions  are  replaced  by  approximations.  In  general,  removal  of  aliasing  errors 
requires  a  decrease  in  the  spatial  step  size  so  that  the  range  of  covered  spatial  frequencies  is  in¬ 
creased.  To  demonstrate  the  effect  of  spatial-frequency  coverage,  a  test  employing  finer  spatial  and 
temporal  sampling  (spatial  step  of  0.0371  mm,  512  x  512  grid  points,  CFL  number  held  at  0.2153) 
was  carried  out  for  the  “bone”  cylinder.  The  required  computation  time  increased  to  136  minutes, 
comparable  to  the  finite-difference  computation  time  of  189  minutes.  Despite  this  increased  com¬ 
putational  cost,  the  time-domain  L2  error  increased  from  0.3072  to  0.4081  for  the  unsmoothed 
case  and  decreased  only  to  0.6841  from  0.8551  for  the  smoothed  case.  This  indicates  that  the 
aliasing  errors  were  not  significantly  reduced  by  this  increase  in  k- space  coverage;  to  achieve  ac¬ 
curacy  comparable  to  finite-difference  computations  for  this  high-contrast  object,  computational 
and  storage  costs  much  greater  than  those  for  the  finite-difference  method  may  be  required. 

Computational  results  for  a  large-scale  two-dimensional  tissue  model  are  shown  in  Fig.  7. 
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Waveforms  computed  by  the  A;-space  (four  points  per  minimum  wavelength,  CFL  number  0.5, 
no  smoothing)  and  the  finite-difference  (ten  points  per  minimum  wavelength,  CFL  number  0.25) 
were  recorded  at  130-element  apertures  composed  of  simulated  point  receivers  separated  by  a 
pitch  of  0.21  mm.  The  results  produced  by  the  finite-difference  method  and  the  /c-space  method 
are  visually  indistinguishable.  However,  despite  the  reduced  grid  size  and  limited  computations 
employed  for  the  finite-difference  method,  the  fc-space  method  was  more  efficient  by  about  a  factor 
of  four:  the  required  CPU  time  for  the  A;-space  method  was  1 .28  hours,  while  the  corresponding 
time  for  the  finite-difference  time-domain  method  was  5.14  hours.  This  discrepancy  in  efficiency 
is  even  more  impressive  when  note  is  made  that  the  A:-space  method  using  4  points  per  minimum 
wavelength  provides  significantly  higher  accuracy  than  the  finite-difference  method  using  10  points 
per  minimum  wavelength  (see  Fig.  3).  Thus,  the  present  A;-space  method  may  be  an  appropriate 
replacement  for  finite-difference  methods  previously  employed  to  compute  propagation  through 
large-scale  soft-tissue  models  [23]— [25]. 
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Y.  Extensions  to  the  A  -Space  Method 


The  present  method  can  be  extended  in  a  number  of  ways  to  increase  its  range  of  applicability  in 
computations  of  ultrasound-tissue  interactions. 

Absorption  effects  could  be  added  to  the  present  algorithm  in  one  of  several  ways.  The  most 
straightforward  method  for  including  absorption  is  to  include  an  ad  hoc  damping  term  proportional 
to  df  /dt  in  Eq.  (3)  [3]— [5].  This  approach  yields  absorption  coefficients  proportional  to  the  square 
of  the  frequency.  Alternatively,  the  sound  speed  c  can  be  replaced  with  a  complex  sound  speed 
including  an  imaginary  part  associated  with  absorption;  this  approach  results  in  absorption  coeffi¬ 
cients  directly  proportional  to  the  frequency.  However,  neither  of  these  approaches  has  a  rigorous 
justification  for  use  in  models  of  ultrasound  propagation  in  biological  tissue. 

A  physically  justifiable  approach  for  inclusion  of  absorption  in  the  present  algorithm  is  to 
consider  absorption  associated  with  multiple  relaxation  processes.  The  theoretical  basis  for  this 
approach  is  presented  in  Ref.  [27];  one  implementation  of  this  method  in  a  finite-difference  time- 
domain  algorithm  is  given  in  Ref.  [28].  Since  multiple  relaxation  processes  can  lead  to  a  variety  of 
frequency-dependent  absorption  characteristics,  this  approach  provides  a  possibility  of  modeling 
realistic  frequency-dependent  attenuation  in  tissue  without  introduction  of  nonphysical  dispersion 
or  violation  of  causality.  Following  the  methods  presented  in  Ref.  [28],  absorption  due  to  multiple 
relaxation  processes  can  be  implemented  in  a  computationally  efficient  form.  Possible  alternatives 
include  the  time-causal  absorption  formulation  of  Ref.  [29]. 

Another  possible  extension  to  the  present  method  is  to  incorporate  the  full  elastic  wave  prop¬ 
agation  equations.  This  extension  would  account  for  shear  wave  propagation,  which  could  sub¬ 
stantially  affect  results  for  propagation  models  including  bone  and  other  calcified  tissue.  By  ap¬ 
plying  methods  similar  to  those  outlined  in  Refs.  [6]  and  [13]  to  the  algorithm  described  above,  a 
full  elastic  /c -space  method  incorporating  Fourier-space  evaluation  of  spatial  derivatives  and  a  k-t 
space  propagator  could  be  implemented. 

Boundary  conditions  of  /c-space  and  pseudospectral  methods  are  inherently  periodic,  so  that 
simple  radiation  boundary  conditions  cannot  be  straightforwardly  implemented.  Instead,  the  tech¬ 
nique  of  perfectly  matched  layers  [30]  may  be  used  in  a  manner  that  reduces  all  waves  to  near 
zero  amplitude  at  the  boundaries  [12, 13, 28].  Extension  of  the  present  method  to  include  perfectly 
matched  layers  is  straightforward. 

Finally,  although  all  implementations  shown  in  the  present  paper  are  two-dimensional,  the 
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present  method  can  be  applied  directly  to  three-dimensional  wave  propagation  without  any  for¬ 
mal  changes.  This  is  possible  because  the  k-t  space  Green’s  function  has  an  identical  form  for 
any  number  of  spatial  dimensions  [2].  To  implement  the  present  methods  for  three-dimensional 
computations,  the  algorithm  outlined  above  is  simply  employed  using  three-dimensional  discrete 
Fourier  transforms. 

Computation  times  for  the  fc-space  method  can  easily  be  reduced  by  parallelization.  The  pri¬ 
mary  computational  burden  of  the  method  is  incurred  in  the  multidimensional  fast  Fourier  trans¬ 
forms  (FFT)  taken  at  each  time  step.  Since  FFT’s  can  be  efficiently  executed  on  parallel  processors 
[31, 32],  the  present  fc-space  method  should  scale  efficiently  to  large  problems  that  require  parallel 
processing. 
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VI.  Conclusions 

A  new  fc-space  method  for  computation  of  ultrasonic  wave  propagation  has  been  presented.  The 
method  efficiently  accounts  for  sound  speed  and  density  variations,  and  can  be  extended  to  include 
realistic  absorption  effects  and  absorbing  boundary  conditions.  Three-dimensional  computations 
can  also  be  performed  without  change  to  the  algorithm  as  derived  here. 

Numerical  results  have  shown  that  the  new  &-space  method  provides  superior  stability  and 
accuracy  over  both  a  similar  leapfrog  pseudospectral  method  and  a  fourth-order  space,  second- 
order  time,  finite-difference  method.  Because  of  this  improved  accuracy,  the  new  method  allows 
larger  spatial  and  time  steps  to  be  employed,  so  that  large-scale  multidimensional  computations 
are  more  feasible  using  the  new  method.  Computations  using  a  realistic  two-dimensional  tissue 
model  support  the  finding  that  the  A;-space  method  provides  high  accuracy  and  low  computational 
cost  for  large-scale  computations. 

The  results  also  indicate  that  care  should  be  taken  when  choosing  and  implementing  a  forward 
solver  for  a  particular  scattering  problem.  For  instance,  in  the  present  A;-space  method,  one  can 
suppress  spurious  waves  by  smoothing  of  contrast  functions;  however,  this  smoothing  can  greatly 
increase  the  time-domain  L2  error  in  some  cases.  Likewise,  the  finite-difference  time-domain 
method  employed  here  is  less  accurate  than  the  /c-space  method  in  most  cases  examined  here,  but 
achieved  higher  accuracy  for  a  test  case  with  a  bone-like  scatterer.  In  general,  the  /c-space  method 
proposed  here  should  be  most  applicable  to  large-scale  scattering  problems  involving  low-contrast 
inhomogeneities  such  as  soft  tissue  structures. 
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Figure  1 :  Time  history  of  total  acoustic  pressure  computed  by  the  fc-space  method  for  a  cylinder  of 
2.0  mm  radius  and  fat-mimicking  acoustic  properties.  The  cylinder  boundary  is  sketched  in  white. 
The  first  panel  shows  the  wavefield  impinging  on  the  cylinder  at  time  t  —  0.96  yus  and  subsequent 
panels  (progressing  from  left  to  right  and  top  to  bottom)  show  the  total  wavefield  at  intervals  of 
0.96  /is.  The  acoustic  pressure  is  plotted  in  all  panels  using  a  bipolar  logarithmic  scale  with  a 
60  dB  dynamic  range. 


Figure  2:  Time-domain  comparison  of  accuracy  for  the  A;-space  and  leapfrog  pseudospectral  meth¬ 
ods  as  a  function  of  CFL  number.  Each  test  used  the  “fat”  cylinder  of  2.0  mm  radius  and  a  spatial 
step  size  of  four  points  per  minimum  wavelength. 
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Figure  3:  Time-domain  comparison  of  accuracy  for  the  /c-space  and  finite-difference  time-domain 
methods  as  a  function  of  the  spatial  step  size  in  points  per  minimum  wavelength  (PPW).  Each  test 
used  the  “fat”  cylinder  of  2.0  mm  radius.  CFL  numbers  were  0.5  for  the  A;-space  method  and  0.25 
for  the  finite-difference  time-domain  method. 
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Figure  4:  Computed  waveforms  for  the  “fat”  cylinder  at  a  radius  of  2.5  mm  for  a  cylinder  of  radius 
2.0  mm  and  a  pulse  center  frequency  of  2.5  MHz.  The  acoustic  pressure  is  shown  on  a  bipolar 
logarithmic  scale  with  60  dB  dynamic  range.  The  horizontal  range  of  each  plot  is  360  degrees, 
covering  the  entire  measurement  circle  starting  with  angle  0  (forward  propagation).  The  vertical 
range  of  each  panel  corresponds  to  a  temporal  duration  of  9.  /is,  with  t  =  0  at  the  top  of  each 
plot,  (a)  Unsmoothed  object;  /c-space  solution  with  four  points  per  minimum  wavelength,  L2  error 
0.0347.  (b)  Smoothed  object;  A;-space  solution  with  four  points  per  minimum  wavelength,  L2  error 
0.0414.  (c)  Finite-difference  solution  with  14  points  per  minimum  wavelength,  L2  error  0.0529. 
(d)  Exact  solution. 
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Figure  5:  Computed  waveforms  at  a  radius  of  12.5  mm  for  a  “fat”  cylinder  of  radius  10.0  mm 
and  a  pulse  center  frequency  of  2.5  MHz.  The  acoustic  pressure  is  shown  in  each  panel  using  a 
bipolar  logarithmic  scale  with  a  60  dB  dynamic  range.  The  horizontal  range  of  each  panel  is  360 
degrees  and  the  vertical  range  is  33  ps.  (a)  Unsmoothed  object;  A;-space  solution,  L2  error  0.1 103. 
(b)  Smoothed  object;  A:-space  solution,  L2  error  0.1080.  (c)  Finite-difference  solution,  L2  error 
0.1763.  (d)  Exact  solution. 
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Figure  6:  Computed  pressure  waveforms  at  a  receiver  radius  of  2.5  mm  for  a  “bone”  cylinder  of 
radius  2.0  mm  and  a  pulse  center  frequency  of  2.5  MHz.  The  format  is  the  same  as  in  Fig.  4.  (a) 
Unsmoothed  object;  /c-space  solution,  L2  error  0.3072.  (b)  Smoothed  object;  fc-space  solution,  L2 
error  0.8551.  (c)  Finite-difference  solution,  L2  error  0.0541.  (d)  Exact  solution. 
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Figure  7:  Comparison  of  /c-space  and  finite-difference  methods  for  a  tissue  cross-sectional  model, 
(a)  Chest  wall  cross  section  (taken  from  Ref.  [24]),  with  black  indicating  connective  tissue,  dark 
gray  indicating  muscle,  and  light  gray  indicating  fat.  The  region  is  33.5  mm  wide  and  17.2  mm 
high,  (b)  Transmitted  waveforms  computed  by  the  /c-space  method  using  four  points  per  minimum 
wavelength  and  a  CFL  number  of  0.5,  shown  on  a  bipolar  linear  gray  scale  with  white  indicating 
maximum  positive  pressure  and  black  indicating  maximum  negative  pressure.  The  horizontal  range 
shown  is  27.3  mm  and  is  shown  to  the  same  scale  as  in  (a).  The  vertical  range  is  3.29  //s.  (c) 
Transmitted  waveforms  computed  by  the  finite-difference  time-domain  method  using  10  points 
per  minimum  wavelength  and  a  CFL  number  of  0.25,  shown  using  the  same  format  as  in  (b). 
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Appendix  D 


Abstracts 

“Simulation  of  Ultrasonic  Pulse  Propagation,  Distortion, 
and  Attenuation  in  the  Human  Chest  Wall” 

Presented  at  the  136th  Meeting  of  the  Acoustical 
Acoustical  Society  of  America 


“Time-Domain  Ultrasound  Diffraction  Tomography” 

Presented  at  Forum  Acusticum  99 
[Joint  Meeting:  137th  Meeting  of  the  Acoustical  Society  of  America, 
2nd  Convention  of  the  European  Acoustics  Association, 
and  25th  German  Acoustics  DAGA  Conference] 


Annual  Report  for  DAMD17-98-1-8141,  July  1999 


4pBB4.  Simulation  of  ultrasonic  propagation,  scattering,  and 
attenuation  in  the  human  chest  wall.  T.  Douglas  Mast  (Appl.  Res. 

Lab.*  Penn  State  Univ.,  University  Park,  PA  16802, 
mast@sabine.acs.psu.edu),  Laura  M.  Hinkelman  (Penn  State  Univ., 
University  Park,  PA  16802),  and  Robert  C.  Waag  (Univ.  of  Rochester, 
Rochester,  NY  14627) 

A  finite-difference  time-domain  model  for  ultrasonic  pulse  propaga¬ 
tion  through  soft  tissue  [Mast  el  al>  J.  Acoust.  Soc.  Am.  102,  1 177-1 190 
(1997)]  has  been  extended  to  incorporate  absorption  effects  as  well  as 
longitudinal-wave  propagation  in  cartilage  and  bone.  This  extended 
method  has  been  used  to  simulate  ultrasonic  propagation  through  anatomi¬ 
cally  detailed  chest  wall  models.  The  inhomogeneous  chest  wall  structure 
is  represented  by  two-dimensional  tissue  maps  determined  by  staining 
chest  wall  cross  sections  to  identify  connective  tissue,  muscle,  fat,  carti¬ 
lage,  and  bone,  scanning  the  stained  cross  sections  at  300  dpi,  and  pro¬ 
cessing  the  scanned  images  to  assign  a  tissue  type  to  each  pixel.  Each 
pixel  of  the  tissue  map  is  then  assigned  a  sound  speed,  density,  and  ab¬ 
sorption  value  determined  from  published  measurements  to  be  representa¬ 
tive  of  the  local  tissue  type.  Computational  results  for  wavefront  distortion 
including  amplitude  fluctuations,  arrival  time  fluctuations,  and  waveform 
distortion  show  qualitative  agreement  with  measurements  performed  on 
the  same  specimens  [Hinkelman  et  al. ,  J.  Acoust.  Soc.  Am.  101,  2365- 
2373  (1997)].  Visualization  of  simulated  tissue-ultrasound  interactions  in 
the  chest  wall  shows  possible  mechanisms  for  image  aberration  in 
echocardiography,  including  effects  due  to  reflection  and  diffraction  from 
rib  structures. 

1844  J.  Acoust.  Soc.  Am.,  Voi.  104,  No.  3,  Pt.  2,  September  1998 


2aPAb6.  Time-domain  ultrasound  diffraction  tomography.  T. 

Douglas  Mast  (Appl.  Res.  Lab.,  Penn  State  Univ.,  University  Park,  PA 
16802,  mast@sabine.acs.psu.edu) 

A  quantitative  ultrasonic  imaging  method  employing  time-domain 
scattering  data  is  presented.  This  method  provides  tomographic  images  of 
inhomogeneous  media  using  scattering  measurements  made  on  a  surface 
surrounding  the  medium  of  interest,  e.g.,  on  a  circle  for  two-dimensional 
problems  or  on  a  sphere  for  three-dimensional  problems.  These  scattering 
data  are  used  to  construct  a  time-domain  analog  of  the  far-field  scattering 
operator.  Images  of  compressibility  variations  are  then  reconstructed  using 
a  coherent  combination  of  the  far-field  scattered  waveforms,  delayed  and 
summed  in  a  manner  that  numerically  focuses  on  the  unknown  medium. 
This  approach  is  closely  related  to  synthetic  aperture  imaging;  however, 
unlike  conventional  synthetic-aperture  methods,  the  present  method  pro¬ 
vides  quantitative  reconstructions  of  compressibility  variations,  analogous 
to  frequency-compunded  filtered  backpropagation  images  weighted  by  the 
spectrum  of  the  incident  wave.  Example  reconstructions,  obtained  using 
synthetic  data  for  two-dimensional  scattering  of  wideband  pulses,  show 
that  the  time-domain  reconstruction  method  can  provide  image  quality 
superior  to  single-frequency  reconstructions  for  objects  of  size  and  con¬ 
trast  relevant  to  medical  imaging  problems  such  as  ultrasonic  mammog¬ 
raphy.  Reconstructions  also  illustrate  the  dependence  of  image  quality  on 
the  number  of  incident- wave  insonifications  and  on  the  range  of  scattering 
angles  available  for  measurements. 


1014 


J.  Acoust.  Soc.  Am.,  Vol.  105,  No.  2,  Pt.  2,  February  1999  Joint  Meeting:  ASA/EAA/DEGA 


Appendix  E 


Curriculum  Vitae 
T.  Douglas  Mast 
Principal  Investigator 


Annual  Report  for  DAMD17-98- 1-8141,  July  1999 


T.  Douglas  Mast 


Research  Associate,  Applied  Research  Laboratory, 

The  Pennsylvania  State  University,  P.O.B.  30,  State  College,  PA  16804 
(814)  863-9998  (tel.),  (814)  863-9918  (fax),  mast@sabine.acs.psu.edu 

Research  and  Teaching  Interests 

Physical  acoustics,  ultrasonic  imaging,  wave  propagation  and  scattering  in  inhomogeneous  media, 
inverse  scattering,  flow/sound  interaction,  bioacoustics,  nondestructive  evaluation. 

Education 

Ph.D.  in  Acoustics,  The  Pennsylvania  State  University,  1993.  Thesis:  Physical  Theory  of  Narrow- 
Band  Sounds  Associated  with  Aneurysms.  Advisor:  Allan  D.  Pierce.  GPA:  3.93/4.0. 

Certificate  (comparable  to  B.A.)  in  Music,  The  Naropa  Institute,  1988. 

B.A.  in  Physics  and  Mathematics,  Goshen  College,  1987.  GPA:  3.89/4.0. 

Present  and  Recent  Employment 

Research  Associate,  Applied  Research  Laboratory,  The  Pennsylvania  State  University,  1997-present. 
Postdoctoral  Scholar,  Applied  Research  Laboratory,  The  Pennsylvania  State  University,  1996-1997. 
Postdoctoral  Fellow,  Ultrasound  Research  Laboratory,  University  of  Rochester,  1993-1996. 
Research  Asst.,  Graduate  Program  in  Acoustics,  Pennsylvania  State  University,  1988-1993. 
Research  Asst.,  Turner  Laboratory  of  Precision  X-Ray  Measurements,  Goshen  College,  1986-1987. 
Teaching  Asst.,  Department  of  Physics,  Goshen  College,  1985-1987. 

Honors,  Awards  and  Society  Affiliations 

Listed  in  Who’s  Who  in  Science  and  Engineering. 

Kenneth  E.  Simowitz  Memorial  Award,  The  Pennsylvania  State  University,  1996. 

F.  V.  Hunt  Fellowship,  Acoustical  Society  of  America,  1994-1995. 

Kenneth  E.  Simowitz  Memorial  Citation,  The  Pennsylvania  State  University,  1992. 

General  Electric  Teaching  Incentive  Loan,  1990. 

Turner  Laboratory  Fellowship,  Goshen  College,  1986. 

Member  of  Acoustical  Society  of  America. 

Member  of  American  Institute  of  Ultrasound  in  Medicine. 

Member  of  Institute  of  Electrical  and  Electronics  Engineers  (Ultrasonics,  Ferroelectrics,  and  Fre¬ 
quency  Control  Society). 

Scientific  Publications 

Mast,  T.  D.  and  Gordon,  G.  A.,  “Quantitative  flaw  reconstruction  from  ultrasonic  surface  wave- 
fields  measured  by  electronic  speckle  pattern  interferometry,”  submitted  to  IEEE  Trans.  Ultrason. 
Ferroelect.  Freq.  Control  (1999). 

Mast,  T.  D.,  “Wideband  quantitative  ultrasonic  imaging  by  time-domain  diffraction  tomography,” 
submitted  to  J.  Acoust.  Soc.  Am.  (1999). 


Mast,  T.  D.,  Hinkelman,  L.  M.,  Metlay,  L.  A.,  Orr,  M.  J.,  and  Waag,  R.  C.,  “Simulation  of 
ultrasonic  pulse  propagation,  distortion,  and  attenuation  in  the  human  chest  wall,”  submitted  to 
J.  Acoust.  Soc.  Am.  (1999). 

Souriau,  L.  P.,  Mast,  T.  D.,  Liu,  D.-L.  D.,  Nachman,  A.  I.,  and  Waag,  R.  C.,  “A  new  k- space 
method  for  large-scale  models  of  wave  propagation  in  tissue,”  submitted  to  IEEE  Trans.  Ultrason. 
Ferroelect.  Freq.  Control  (1999). 

Mast,  T.  D.,  Swanson,  D.  C.,  Mahon,  M.  P.  and  Norris,  D.  E.,  “Resolution  of  multipath  outdoor 
sound  propagation  using  spread  spectrum  signals,”  submitted  to  J.  Acoust.  Soc.  Am.  (1998). 

Gordon,  G.  A.  and  Mast,  T.  D.,  “Wide-area  imaging  of  ultrasonic  Lamb  wave  fields  by  electronic 
speckle  pattern  interferometry,”  Proc.  SPIE  3586,  297-309  (1999). 

Myers,  L.  F,  Lovette,  M.,  Kilgus,  C.  C.,  Giannini,  J.  A.,  Swanson,  D.  C.,  Reichard,  K.  M.,  Ma¬ 
hon,  M.  P.,  and  Mast,  T.  D.,  “Java-based  information  system  for  wayside  sensing  and  control,” 
Proceedings  of  the  IEEE/ ASME  Joint  Railroad  Conference ,  135-147  (1998). 

Hinkelman,  L.  M.,  Mast,  T.  D.,  Metlay,  L.  A.,  and  Waag,  R.  C.,  “The  effect  of  abdominal  wall 
morphology  on  ultrasonic  pulse  distortion.  Part  I:  Measurements,”  J.  Acoust.  Soc.  Am.  104,  3635- 
3649  (1998). 

Mast,  T.  D.,  Hinkelman,  L.  M.,  Orr,  M.  J.,  and  Waag,  R.  C.,  “The  effect  of  abdominal  wall 
morphology  on  ultrasonic  pulse  distortion.  Part  II:  Simulations,”  J.  Acoust.  Soc.  Am.  104,  3650- 
3664  (1998). 

Mast,  T.  D.,  Hinkelman,  L.  M.,  Orr,  M.  J.,  Sparrow,  V.  W.,  and  Waag,  R.  C.,  Erratum:  “Simulation 
of  ultrasonic  pulse  propagation  through  abdominal  wall,”  [J.  Acoust.  Soc.  Am.  102,  1177-1190 
(1997)],  J.  Acoust.  Soc.  Am.  104,  1124-1125  (1998). 

Jansson,  T.  T.,  Mast  T.  D.,  and  Waag,  R.  C.,  “Measurements  of  differential  scattering  cross-section 
using  a  ring  transducer,”  J.  Acoust.  Soc.  Am.  103,  3169-3179  (1998). 

Mast,  T.  D.,  Nachman,  A.  I.,  Liu,  D.-L.,  and  Waag,  R.  C.,  “Quantitative  imaging  with  eigen¬ 
functions  of  the  scattering  operator,”  1997  IEEE  Ultrasonics  Symposium  Proceedings,  Vol.  2,  pp. 
1507-1510. 

Hinkelman,  L.  M.,  Mast,  T.  D.,  Orr,  M.  J.,  and  Waag,  R.  C.,  “Effects  of  abdominal  wall  morphology 
on  ultrasonic  pulses,”  1997  IEEE  Ultrasonics  Symposium  Proceedings,  Vol.  2,  pp.  1493-1496. 

Mast,  T.  D.,  Nachman,  A.  I.,  and  Waag,  R.  C.,  “Focusing  and  imaging  using  eigenfunctions  of  the 
scattering  operator,”  J.  Acoust.  Soc.  Am.  102,  715-725  (1997). 

Mast,  T.  D.,  Hinkelman,  L.  M.,  Orr,  M.  J.,  Sparrow,  V.  W.,  and  Waag,  R.  C.,  “Simulation  of 
ultrasonic  pulse  propagation  through  abdominal  wall,”  J.  Acoust.  Soc.  Am.  102,  1177-1190  (1997). 

Mast,  T.  D.  and  Waag,  R.  C.,  “Wave  space  resolution  in  ultrasonic  scattering  measurements,”  J. 
Acoust.  Soc.  Am.  98,  3050-3058  (1995). 

Mast,  T.  D.  and  Pierce,  A.  D.,  “A  theory  of  aneurysm  sounds,”  J.  Biomech.  28,  1045-1053  (1995). 

Mast,  T.  D.  and  Pierce,  A.  D.,  “Describing-function  theory  for  flow  excitation  of  resonators,”  J. 
Acoust.  Soc.  Am.  97,  163-172  (1995). 

Mast,  T.  D.,  “Limit  cycles  of  flow-excited  resonators:  a  describing- function  analysis,”  in  Struc¬ 
tural  Acoustics,  Scattering,  and  Propagation:  Theoretical  and  Computational  Acoustics — Volume 
I,  Edited  by  J.  E.  Ffowcs  Williams,  D.  Lee,  and  A.  D.  Pierce,  (River  Edge,  New  Jersey:  World 
Scientific,  1994),  pp.  389-403. 

Mast,  T.  D.  and  Pierce,  A.  D.,  “Flow-induced  sounds  associated  with  aneurysms,”  in  Flow  Noise 
Modeling,  Measurement,  and  Control,  edited  by  T.  M.  Farabee,  W.  L.  Keith,  and  R.  M.  Lueptow 
(New  York:  American  Society  of  Mechanical  Engineers,  1991),  pp.  129-134. 


